|
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 ! ***************************************************************************** 00015 MODULE md_vel_utils 00016 USE atomic_kind_list_types, ONLY: atomic_kind_list_type 00017 USE atomic_kind_types, ONLY: atomic_kind_type,& 00018 get_atomic_kind,& 00019 get_atomic_kind_set 00020 USE cell_types, ONLY: cell_type 00021 USE cp_linked_list_val, ONLY: cp_sll_val_next,& 00022 cp_sll_val_type 00023 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 00024 cp_print_key_unit_nr 00025 USE cp_para_types, ONLY: cp_para_env_type 00026 USE cp_subsys_types, ONLY: cp_subsys_get,& 00027 cp_subsys_type 00028 USE cp_units, ONLY: cp_unit_from_cp2k 00029 USE extended_system_types, ONLY: npt_info_type 00030 USE f77_blas 00031 USE force_env_methods, ONLY: force_env_calc_energy_force 00032 USE force_env_types, ONLY: force_env_get,& 00033 force_env_type 00034 USE force_env_utils, ONLY: force_env_rattle,& 00035 force_env_shake 00036 USE global_types, ONLY: global_environment_type 00037 USE input_constants, ONLY: & 00038 npe_f_ensemble, npe_i_ensemble, nph_uniaxial_damped_ensemble, & 00039 nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, & 00040 reftraj_ensemble, use_perd_none, use_perd_x, use_perd_xy, & 00041 use_perd_xyz, use_perd_xz, use_perd_y, use_perd_yz, use_perd_z 00042 USE input_cp2k_binary_restarts, ONLY: read_binary_velocities 00043 USE input_cp2k_restarts, ONLY: update_subsys 00044 USE input_section_types, ONLY: section_vals_get,& 00045 section_vals_get_subs_vals,& 00046 section_vals_list_get,& 00047 section_vals_type,& 00048 section_vals_val_get 00049 USE input_val_types, ONLY: val_get,& 00050 val_type 00051 USE kinds, ONLY: default_string_length,& 00052 dp 00053 USE mathlib, ONLY: diamat_all 00054 USE md_ener_types, ONLY: md_ener_type 00055 USE md_environment_types, ONLY: get_md_env,& 00056 md_environment_type 00057 USE mol_kind_new_list_types, ONLY: mol_kind_new_list_type 00058 USE molecule_kind_types, ONLY: fixd_constraint_type,& 00059 get_molecule_kind,& 00060 get_molecule_kind_set,& 00061 molecule_kind_type 00062 USE parallel_rng_types, ONLY: next_random_number 00063 USE particle_list_types, ONLY: particle_list_type 00064 USE particle_types, ONLY: particle_type 00065 USE physcon, ONLY: kelvin 00066 USE shell_opt, ONLY: optimize_shell_core 00067 USE shell_potential_types, ONLY: shell_kind_type 00068 USE simpar_types, ONLY: simpar_type 00069 USE thermal_region_types, ONLY: thermal_region_type,& 00070 thermal_regions_type 00071 USE timings, ONLY: timeset,& 00072 timestop 00073 #include "cp_common_uses.h" 00074 00075 IMPLICIT NONE 00076 00077 PRIVATE 00078 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'md_vel_utils' 00079 00080 00081 PUBLIC :: temperature_control,& 00082 comvel_control,& 00083 angvel_control,& 00084 setup_velocities 00085 00086 CONTAINS 00087 00088 ! ***************************************************************************** 00095 SUBROUTINE compute_rcom(part,is_fixed,rcom) 00096 TYPE(particle_type), DIMENSION(:), 00097 POINTER :: part 00098 INTEGER, DIMENSION(:), INTENT(IN) :: is_fixed 00099 REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: rcom 00100 00101 CHARACTER(len=*), PARAMETER :: routineN = 'compute_rcom', 00102 routineP = moduleN//':'//routineN 00103 00104 INTEGER :: i 00105 REAL(KIND=dp) :: denom, mass 00106 TYPE(atomic_kind_type), POINTER :: atomic_kind 00107 00108 rcom(:) = 0.0_dp 00109 denom = 0.0_dp 00110 DO i = 1, SIZE(part) 00111 atomic_kind => part(i)%atomic_kind 00112 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 00113 SELECT CASE(is_fixed(i)) 00114 CASE(use_perd_x,use_perd_y,use_perd_z,use_perd_xy,use_perd_xz,use_perd_yz,use_perd_none) 00115 rcom(1) = rcom(1) + part(i)%r(1) * mass 00116 rcom(2) = rcom(2) + part(i)%r(2) * mass 00117 rcom(3) = rcom(3) + part(i)%r(3) * mass 00118 denom = denom + mass 00119 END SELECT 00120 END DO 00121 rcom = rcom/denom 00122 00123 END SUBROUTINE compute_rcom 00124 00125 ! ***************************************************************************** 00132 SUBROUTINE compute_vcom(part,is_fixed,vcom,ecom) 00133 TYPE(particle_type), DIMENSION(:), 00134 POINTER :: part 00135 INTEGER, DIMENSION(:), INTENT(IN) :: is_fixed 00136 REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: vcom 00137 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: ecom 00138 00139 CHARACTER(len=*), PARAMETER :: routineN = 'compute_vcom', 00140 routineP = moduleN//':'//routineN 00141 00142 INTEGER :: i 00143 REAL(KIND=dp) :: denom, mass 00144 TYPE(atomic_kind_type), POINTER :: atomic_kind 00145 00146 vcom = 0.0_dp 00147 denom = 0.0_dp 00148 DO i = 1, SIZE(part) 00149 atomic_kind => part(i)%atomic_kind 00150 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 00151 IF (mass.NE.0.0) THEN 00152 SELECT CASE(is_fixed(i)) 00153 CASE(use_perd_x,use_perd_y,use_perd_z,use_perd_xy,use_perd_xz,use_perd_yz,use_perd_none) 00154 vcom(1) = vcom(1) + part(i)%v(1) * mass 00155 vcom(2) = vcom(2) + part(i)%v(2) * mass 00156 vcom(3) = vcom(3) + part(i)%v(3) * mass 00157 denom = denom + mass 00158 END SELECT 00159 END IF 00160 END DO 00161 vcom = vcom/denom 00162 IF (PRESENT(ecom)) THEN 00163 ecom = 0.5_dp*denom*SUM(vcom*vcom) 00164 END IF 00165 00166 END SUBROUTINE compute_vcom 00167 00168 ! ***************************************************************************** 00175 SUBROUTINE clone_core_shell_vel(part,shell_part,core_part) 00176 TYPE(particle_type), DIMENSION(:), 00177 POINTER :: part, shell_part, core_part 00178 00179 CHARACTER(len=*), PARAMETER :: routineN = 'clone_core_shell_vel', 00180 routineP = moduleN//':'//routineN 00181 00182 INTEGER :: i 00183 LOGICAL :: is_shell 00184 TYPE(atomic_kind_type), POINTER :: atomic_kind 00185 00186 DO i = 1, SIZE(part) 00187 atomic_kind => part(i)%atomic_kind 00188 CALL get_atomic_kind(atomic_kind=atomic_kind, shell_active=is_shell) 00189 IF (is_shell) THEN 00190 shell_part( part(i)%shell_index )%v(:) = part(i)%v(:) 00191 core_part( part(i)%shell_index )%v(:) = part(i)%v(:) 00192 END IF 00193 END DO 00194 00195 END SUBROUTINE clone_core_shell_vel 00196 00197 ! ***************************************************************************** 00205 FUNCTION compute_ekin(part,ireg) RESULT(ekin) 00206 TYPE(particle_type), DIMENSION(:), 00207 POINTER :: part 00208 INTEGER, INTENT(IN), OPTIONAL :: ireg 00209 REAL(KIND=dp) :: ekin 00210 00211 CHARACTER(len=*), PARAMETER :: routineN = 'compute_ekin', 00212 routineP = moduleN//':'//routineN 00213 00214 INTEGER :: i 00215 REAL(KIND=dp) :: mass 00216 TYPE(atomic_kind_type), POINTER :: atomic_kind 00217 00218 NULLIFY(atomic_kind) 00219 ekin = 0.0_dp 00220 IF(PRESENT(ireg)) THEN 00221 DO i = 1, SIZE(part) 00222 IF(part(i)%t_region_index==ireg) THEN 00223 atomic_kind => part(i)%atomic_kind 00224 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 00225 ekin = ekin + 0.5_dp * mass * SUM(part(i)%v(:) * part(i)%v(:)) 00226 END IF 00227 END DO 00228 ELSE 00229 DO i = 1, SIZE(part) 00230 atomic_kind => part(i)%atomic_kind 00231 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 00232 ekin = ekin + 0.5_dp * mass * SUM(part(i)%v(:) * part(i)%v(:)) 00233 END DO 00234 END IF 00235 00236 END FUNCTION compute_ekin 00237 00238 ! ***************************************************************************** 00246 SUBROUTINE rescale_vel(part,simpar,ekin,vcom,ireg,nfree,temp) 00247 TYPE(particle_type), DIMENSION(:), 00248 POINTER :: part 00249 TYPE(simpar_type), POINTER :: simpar 00250 REAL(KIND=dp), INTENT(INOUT) :: ekin 00251 REAL(KIND=dp), DIMENSION(3), 00252 INTENT(INOUT), OPTIONAL :: vcom 00253 INTEGER, INTENT(IN), OPTIONAL :: ireg, nfree 00254 REAL(KIND=dp), INTENT(IN), OPTIONAL :: temp 00255 00256 CHARACTER(len=*), PARAMETER :: routineN = 'rescale_vel', 00257 routineP = moduleN//':'//routineN 00258 00259 INTEGER :: i, my_ireg, my_nfree 00260 REAL(KIND=dp) :: factor, my_temp 00261 00262 IF(PRESENT(ireg).AND.PRESENT(nfree).AND.PRESENT(temp)) THEN 00263 my_ireg = ireg 00264 my_nfree = nfree 00265 my_temp = temp 00266 ELSEIF(PRESENT(nfree)) THEN 00267 my_ireg = 0 00268 my_nfree = nfree 00269 my_temp = simpar%temp_ext 00270 ELSE 00271 my_ireg = 0 00272 my_nfree = simpar%nfree 00273 my_temp = simpar%temp_ext 00274 END IF 00275 IF (my_nfree/=0) THEN 00276 factor = my_temp / ( 2.0_dp * ekin ) * REAL(my_nfree,KIND=dp) 00277 ELSE 00278 factor = 0.0_dp 00279 ENDIF 00280 ! Note: 00281 ! this rescaling is still wrong, it should take the masses into account 00282 ! rescaling is generally not correct, so needs fixing 00283 ekin = ekin * factor 00284 factor = SQRT(factor) 00285 IF(PRESENT(ireg)) THEN 00286 DO i = 1, SIZE(part) 00287 IF(part(i)%t_region_index == my_ireg) part(i)%v(:) = factor*part(i)%v(:) 00288 END DO 00289 ELSE 00290 DO i = 1, SIZE(part) 00291 part(i)%v(:) = factor*part(i)%v(:) 00292 END DO 00293 IF (PRESENT(vcom)) THEN 00294 vcom = factor*vcom 00295 END IF 00296 END IF 00297 00298 END SUBROUTINE rescale_vel 00299 00300 ! ***************************************************************************** 00306 SUBROUTINE rescale_vel_region(part,md_env,simpar,error) 00307 00308 TYPE(particle_type), DIMENSION(:), 00309 POINTER :: part 00310 TYPE(md_environment_type), POINTER :: md_env 00311 TYPE(simpar_type), POINTER :: simpar 00312 TYPE(cp_error_type), INTENT(inout) :: error 00313 00314 CHARACTER(LEN=*), PARAMETER :: routineN = 'rescale_vel_region', 00315 routineP = moduleN//':'//routineN 00316 00317 INTEGER :: ireg, nfree, nfree0, 00318 nfree_done 00319 REAL(KIND=dp) :: ekin, temp 00320 TYPE(thermal_region_type), POINTER :: t_region 00321 TYPE(thermal_regions_type), POINTER :: thermal_regions 00322 00323 NULLIFY(thermal_regions, t_region) 00324 00325 CALL get_md_env(md_env,thermal_regions=thermal_regions,error=error) 00326 nfree_done = 0 00327 DO ireg = 1,thermal_regions%nregions 00328 NULLIFY(t_region) 00329 t_region => thermal_regions%thermal_region(ireg) 00330 nfree = t_region%npart*3 00331 ekin = compute_ekin(part,ireg) 00332 temp = t_region%temp_expected 00333 CALL rescale_vel(part,simpar,ekin,ireg=ireg,nfree=nfree,temp=temp) 00334 nfree_done = nfree_done+nfree 00335 ekin = compute_ekin(part,ireg) 00336 temp = 2.0_dp* ekin/REAL(nfree,dp)*kelvin 00337 t_region%temperature = temp 00338 END DO 00339 nfree0 = simpar%nfree - nfree_done 00340 IF(nfree0>0) THEN 00341 ekin = compute_ekin(part,0) 00342 CALL rescale_vel(part,simpar,ekin,ireg=0,nfree=nfree0,temp=simpar%temp_ext) 00343 ekin = compute_ekin(part,0) 00344 temp = 2.0_dp* ekin/REAL(nfree0,dp)*kelvin 00345 thermal_regions%temp_reg0 = temp 00346 END IF 00347 END SUBROUTINE rescale_vel_region 00348 00349 ! ***************************************************************************** 00356 SUBROUTINE subtract_vcom(part,is_fixed,vcom) 00357 TYPE(particle_type), DIMENSION(:), 00358 POINTER :: part 00359 INTEGER, DIMENSION(:), INTENT(IN) :: is_fixed 00360 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: vcom 00361 00362 CHARACTER(len=*), PARAMETER :: routineN = 'subtract_vcom', 00363 routineP = moduleN//':'//routineN 00364 00365 INTEGER :: i 00366 00367 DO i = 1, SIZE(part) 00368 SELECT CASE(is_fixed(i)) 00369 CASE(use_perd_x) 00370 part(i)%v(2) = part(i)%v(2) - vcom(2) 00371 part(i)%v(3) = part(i)%v(3) - vcom(3) 00372 CASE(use_perd_y) 00373 part(i)%v(1) = part(i)%v(1) - vcom(1) 00374 part(i)%v(3) = part(i)%v(3) - vcom(3) 00375 CASE(use_perd_z) 00376 part(i)%v(1) = part(i)%v(1) - vcom(1) 00377 part(i)%v(2) = part(i)%v(2) - vcom(2) 00378 CASE(use_perd_xy) 00379 part(i)%v(3) = part(i)%v(3) - vcom(3) 00380 CASE(use_perd_xz) 00381 part(i)%v(2) = part(i)%v(2) - vcom(2) 00382 CASE(use_perd_yz) 00383 part(i)%v(1) = part(i)%v(1) - vcom(1) 00384 CASE(use_perd_none) 00385 part(i)%v(:) = part(i)%v(:) - vcom(:) 00386 END SELECT 00387 END DO 00388 END SUBROUTINE subtract_vcom 00389 00390 ! ***************************************************************************** 00397 SUBROUTINE compute_vang(part,is_fixed,rcom,vang,error) 00398 TYPE(particle_type), DIMENSION(:), 00399 POINTER :: part 00400 INTEGER, DIMENSION(:), INTENT(IN) :: is_fixed 00401 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rcom 00402 REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: vang 00403 TYPE(cp_error_type), INTENT(inout) :: error 00404 00405 CHARACTER(len=*), PARAMETER :: routineN = 'compute_vang', 00406 routineP = moduleN//':'//routineN 00407 00408 INTEGER :: i 00409 REAL(KIND=dp) :: mass, proj 00410 REAL(KIND=dp), DIMENSION(3) :: evals, mang, r 00411 REAL(KIND=dp), DIMENSION(3, 3) :: iner 00412 TYPE(atomic_kind_type), POINTER :: atomic_kind 00413 00414 NULLIFY(atomic_kind) 00415 mang(:) = 0.0_dp 00416 iner(:,:) = 0.0_dp 00417 DO i=1,SIZE(part) 00418 ! compute angular momentum and inertia tensor 00419 SELECT CASE(is_fixed(i)) 00420 CASE(use_perd_x,use_perd_y,use_perd_z,use_perd_xy,use_perd_xz,use_perd_yz,use_perd_none) 00421 r(:) = part(i)%r(:) - rcom(:) 00422 atomic_kind => part(i)%atomic_kind 00423 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 00424 mang(1) = mang(1) + mass*(r(2)*part(i)%v(3) - r(3)*part(i)%v(2)) 00425 mang(2) = mang(2) + mass*(r(3)*part(i)%v(1) - r(1)*part(i)%v(3)) 00426 mang(3) = mang(3) + mass*(r(1)*part(i)%v(2) - r(2)*part(i)%v(1)) 00427 00428 iner(1,1) = iner(1,1) + mass*(r(2)*r(2) + r(3)*r(3)) 00429 iner(2,2) = iner(2,2) + mass*(r(3)*r(3) + r(1)*r(1)) 00430 iner(3,3) = iner(3,3) + mass*(r(1)*r(1) + r(2)*r(2)) 00431 00432 iner(1,2) = iner(1,2) - mass*r(1)*r(2) 00433 iner(2,3) = iner(2,3) - mass*r(2)*r(3) 00434 iner(3,1) = iner(3,1) - mass*r(3)*r(1) 00435 END SELECT 00436 END DO 00437 iner(2,1) = iner(1,2) 00438 iner(3,2) = iner(2,3) 00439 iner(1,3) = iner(3,1) 00440 00441 ! Take the safest route, i.e. diagonalize the inertia tensor and solve 00442 ! the angular velocity only with the non-zero eigenvalues. A plain inversion 00443 ! would fail for linear molecules. 00444 CALL diamat_all(iner, evals,error=error) 00445 00446 vang(:) = 0.0_dp 00447 DO i=1,3 00448 IF (evals(i) > 0.0_dp) THEN 00449 proj = SUM(iner(:,i)*mang)/evals(i) 00450 vang(1) = vang(1) + proj*iner(1,i) 00451 vang(2) = vang(2) + proj*iner(2,i) 00452 vang(3) = vang(3) + proj*iner(3,i) 00453 END IF 00454 END DO 00455 00456 END SUBROUTINE compute_vang 00457 00458 ! ***************************************************************************** 00465 SUBROUTINE subtract_vang(part,is_fixed,rcom,vang) 00466 TYPE(particle_type), DIMENSION(:), 00467 POINTER :: part 00468 INTEGER, DIMENSION(:), INTENT(IN) :: is_fixed 00469 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rcom, vang 00470 00471 CHARACTER(len=*), PARAMETER :: routineN = 'subtract_vang', 00472 routineP = moduleN//':'//routineN 00473 00474 INTEGER :: i 00475 REAL(KIND=dp), DIMENSION(3) :: r 00476 00477 DO i=1,SIZE(part) 00478 r(:) = part(i)%r(:) - rcom(:) 00479 SELECT CASE(is_fixed(i)) 00480 CASE(use_perd_x) 00481 part(i)%v(2) = part(i)%v(2) - (vang(3)*r(1) - vang(1)*r(3)) 00482 part(i)%v(3) = part(i)%v(3) - (vang(1)*r(2) - vang(2)*r(1)) 00483 CASE(use_perd_y) 00484 part(i)%v(1) = part(i)%v(1) - (vang(2)*r(3) - vang(3)*r(2)) 00485 part(i)%v(3) = part(i)%v(3) - (vang(1)*r(2) - vang(2)*r(1)) 00486 CASE(use_perd_z) 00487 part(i)%v(1) = part(i)%v(1) - (vang(2)*r(3) - vang(3)*r(2)) 00488 part(i)%v(2) = part(i)%v(2) - (vang(3)*r(1) - vang(1)*r(3)) 00489 CASE(use_perd_xy) 00490 part(i)%v(3) = part(i)%v(3) - (vang(1)*r(2) - vang(2)*r(1)) 00491 CASE(use_perd_xz) 00492 part(i)%v(2) = part(i)%v(2) - (vang(3)*r(1) - vang(1)*r(3)) 00493 CASE(use_perd_yz) 00494 part(i)%v(1) = part(i)%v(1) - (vang(2)*r(3) - vang(3)*r(2)) 00495 CASE(use_perd_none) 00496 part(i)%v(1) = part(i)%v(1) - (vang(2)*r(3) - vang(3)*r(2)) 00497 part(i)%v(2) = part(i)%v(2) - (vang(3)*r(1) - vang(1)*r(3)) 00498 part(i)%v(3) = part(i)%v(3) - (vang(1)*r(2) - vang(2)*r(1)) 00499 END SELECT 00500 END DO 00501 00502 END SUBROUTINE subtract_vang 00503 00504 ! ***************************************************************************** 00512 SUBROUTINE initialize_velocities(simpar, part, force_env, globenv, md_env, & 00513 molecule_kinds, label, print_section, subsys_section, shell_present, shell_part, & 00514 core_part, force_rescaling, para_env, error) 00515 00516 TYPE(simpar_type), POINTER :: simpar 00517 TYPE(particle_type), DIMENSION(:), 00518 POINTER :: part 00519 TYPE(force_env_type), POINTER :: force_env 00520 TYPE(global_environment_type), POINTER :: globenv 00521 TYPE(md_environment_type), POINTER :: md_env 00522 TYPE(mol_kind_new_list_type), POINTER :: molecule_kinds 00523 CHARACTER(LEN=*), INTENT(IN) :: label 00524 TYPE(section_vals_type), POINTER :: print_section, subsys_section 00525 LOGICAL, INTENT(IN) :: shell_present 00526 TYPE(particle_type), DIMENSION(:), 00527 POINTER :: shell_part, core_part 00528 LOGICAL, INTENT(IN) :: force_rescaling 00529 TYPE(cp_para_env_type), POINTER :: para_env 00530 TYPE(cp_error_type), INTENT(INOUT) :: error 00531 00532 CHARACTER(LEN=*), PARAMETER :: routineN = 'initialize_velocities', 00533 routineP = moduleN//':'//routineN 00534 00535 CHARACTER(LEN=default_string_length) :: my_format 00536 INTEGER :: handle, i, ifixd, 00537 imolecule_kind, iw, natoms, 00538 nshell, num_x, shell_index, 00539 stat 00540 INTEGER, ALLOCATABLE, DIMENSION(:) :: is_fixed 00541 LOGICAL :: atomvel_explicit, atomvel_read, corevel_explicit, 00542 corevel_read, failure, is_ok, rescale_regions, shellvel_explicit, 00543 shellvel_read 00544 REAL(KIND=dp) :: ecom, ekin, fac_massc, 00545 fac_masss, mass, mass_tot, 00546 temp, tmp_r1 00547 REAL(KIND=dp), DIMENSION(3) :: rcom, vang, vc, vcom, vs 00548 REAL(KIND=dp), DIMENSION(:), POINTER :: vel 00549 TYPE(atomic_kind_type), POINTER :: atomic_kind 00550 TYPE(cell_type), POINTER :: cell 00551 TYPE(cp_logger_type), POINTER :: logger 00552 TYPE(cp_sll_val_type), POINTER :: atom_list, core_list, 00553 shell_list 00554 TYPE(fixd_constraint_type), 00555 DIMENSION(:), POINTER :: fixd_list 00556 TYPE(molecule_kind_type), DIMENSION(:), 00557 POINTER :: molecule_kind_set 00558 TYPE(molecule_kind_type), POINTER :: molecule_kind 00559 TYPE(section_vals_type), POINTER :: atomvel_section, 00560 corevel_section, 00561 shellvel_section 00562 TYPE(shell_kind_type), POINTER :: shell 00563 TYPE(thermal_regions_type), POINTER :: thermal_regions 00564 TYPE(val_type), POINTER :: val 00565 00566 CALL timeset(routineN,handle) 00567 00568 ! Initializing parameters 00569 failure = .FALSE. 00570 NULLIFY (atomic_kind, fixd_list, atom_list, shell_list, core_list, logger, molecule_kind) 00571 NULLIFY (molecule_kind_set, shell, val) 00572 NULLIFY (atomvel_section,shellvel_section, corevel_section) 00573 atomvel_read = .FALSE. 00574 corevel_read = .FALSE. 00575 shellvel_read = .FALSE. 00576 00577 atomvel_section => section_vals_get_subs_vals(subsys_section,"VELOCITY",error=error) 00578 shellvel_section => section_vals_get_subs_vals(subsys_section,"SHELL_VELOCITY",error=error) 00579 corevel_section => section_vals_get_subs_vals(subsys_section,"CORE_VELOCITY",error=error) 00580 00581 ! Logging 00582 logger => cp_error_get_logger(error) 00583 iw=cp_print_key_unit_nr(logger,print_section,"PROGRAM_RUN_INFO",extension=".log",error=error) 00584 IF (iw>0) THEN 00585 num_x = (79-LEN_TRIM(ADJUSTL(label))-2)/2 00586 WRITE(my_format,'(A,I0,A,I0,A)')'(1X,',num_x,'("*"),1X,A,1X,',num_x,'("*"))' 00587 WRITE ( iw, TRIM(my_format))TRIM(ADJUSTL(label)) 00588 END IF 00589 00590 ! Initializing parameters 00591 natoms = SIZE(part) 00592 nshell = 0 00593 00594 ! Core-Shell Model 00595 IF (shell_present) THEN 00596 CPPostcondition(ASSOCIATED(core_part),cp_failure_level,routineP,error,failure) 00597 CPPostcondition(ASSOCIATED(shell_part),cp_failure_level,routineP,error,failure) 00598 nshell = SIZE(shell_part) 00599 END IF 00600 00601 ! Build a list of all fixed atoms (if any) 00602 ALLOCATE (is_fixed(natoms),STAT=stat) 00603 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00604 00605 is_fixed = use_perd_none 00606 molecule_kind_set => molecule_kinds%els 00607 DO imolecule_kind=1,molecule_kinds%n_els 00608 molecule_kind => molecule_kind_set(imolecule_kind) 00609 CALL get_molecule_kind(molecule_kind=molecule_kind,fixd_list=fixd_list) 00610 IF (ASSOCIATED(fixd_list)) THEN 00611 DO ifixd=1,SIZE(fixd_list) 00612 IF (.NOT.fixd_list(ifixd)%restraint%active) is_fixed(fixd_list(ifixd)%fixd) = fixd_list(ifixd)%itype 00613 END DO 00614 END IF 00615 END DO 00616 00617 ! Compute the total mass when needed 00618 IF ( simpar%ensemble == nph_uniaxial_ensemble .OR.& 00619 simpar%ensemble == nph_uniaxial_damped_ensemble ) THEN 00620 mass_tot = 0.0_dp 00621 DO i = 1, natoms 00622 atomic_kind => part(i)%atomic_kind 00623 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 00624 mass_tot = mass_tot + mass 00625 END DO 00626 simpar%v_shock = simpar%v_shock * SQRT ( mass_tot ) 00627 END IF 00628 00629 ! Read velocities from binary restart file if available 00630 CALL read_binary_velocities("",part,force_env%root_section,para_env,& 00631 subsys_section,atomvel_read,error) 00632 CALL read_binary_velocities("SHELL",shell_part,force_env%root_section,para_env,& 00633 subsys_section,shellvel_read,error) 00634 CALL read_binary_velocities("CORE",core_part,force_env%root_section,para_env,& 00635 subsys_section,corevel_read,error) 00636 00637 ! Read or initialize the particle velocities 00638 CALL section_vals_get(atomvel_section,explicit=atomvel_explicit,error=error) 00639 CALL section_vals_get(shellvel_section,explicit=shellvel_explicit,error=error) 00640 CALL section_vals_get(corevel_section,explicit=corevel_explicit,error=error) 00641 CPPostcondition(shellvel_explicit.EQV.corevel_explicit,cp_failure_level,routineP,error,failure) 00642 00643 IF (atomvel_explicit.OR.atomvel_read) THEN 00644 IF (.NOT.atomvel_read) THEN 00645 ! Read the atom velocities if explicitly given in the input file 00646 CALL section_vals_list_get(atomvel_section,"_DEFAULT_KEYWORD_",list=atom_list,error=error) 00647 DO i = 1, natoms 00648 is_ok = cp_sll_val_next(atom_list,val,error=error) 00649 CALL val_get(val,r_vals=vel,error=error) 00650 part(i)%v = vel 00651 END DO 00652 END IF 00653 DO i = 1, natoms 00654 SELECT CASE(is_fixed(i)) 00655 CASE (use_perd_x) 00656 part(i)%v(1) = 0.0_dp 00657 CASE (use_perd_y) 00658 part(i)%v(2) = 0.0_dp 00659 CASE (use_perd_z) 00660 part(i)%v(3) = 0.0_dp 00661 CASE (use_perd_xy) 00662 part(i)%v(1) = 0.0_dp 00663 part(i)%v(2) = 0.0_dp 00664 CASE (use_perd_xz) 00665 part(i)%v(1) = 0.0_dp 00666 part(i)%v(3) = 0.0_dp 00667 CASE (use_perd_yz) 00668 part(i)%v(2) = 0.0_dp 00669 part(i)%v(3) = 0.0_dp 00670 CASE (use_perd_xyz) 00671 part(i)%v = 0.0_dp 00672 END SELECT 00673 END DO 00674 IF (shell_present) THEN 00675 IF (shellvel_explicit) THEN 00676 ! If the atoms positions are given (?) and core and shell velocities are 00677 ! present in the input, read the latter. 00678 CALL section_vals_list_get(shellvel_section,"_DEFAULT_KEYWORD_",list=shell_list,error=error) 00679 CALL section_vals_list_get(corevel_section,"_DEFAULT_KEYWORD_",list=core_list,error=error) 00680 DO i=1,nshell 00681 is_ok = cp_sll_val_next(shell_list,val,error=error) 00682 CALL val_get(val,r_vals=vel,error=error) 00683 shell_part(i)%v = vel 00684 is_ok = cp_sll_val_next(core_list,val,error=error) 00685 CALL val_get(val,r_vals=vel,error=error) 00686 core_part(i)%v = vel 00687 END DO 00688 ELSE 00689 IF (.NOT.(shellvel_read.AND.corevel_read)) THEN 00690 ! Otherwise, just copy atom velocties into shell and core velocities. 00691 CALL clone_core_shell_vel(part,shell_part,core_part) 00692 END IF 00693 END IF 00694 END IF 00695 00696 ! compute vcom, ecom and ekin 00697 CALL compute_vcom(part,is_fixed,vcom,ecom) 00698 ekin = compute_ekin(part) - ecom 00699 00700 IF(simpar%do_thermal_region) THEN 00701 CALL get_md_env (md_env, thermal_regions=thermal_regions, error=error) 00702 IF(ASSOCIATED(thermal_regions)) THEN 00703 rescale_regions = thermal_regions%force_rescaling 00704 END IF 00705 ELSE 00706 rescale_regions = .FALSE. 00707 END IF 00708 IF (simpar%nfree /= 0 .AND. (force_rescaling .OR. rescale_regions)) THEN 00709 IF(simpar%do_thermal_region) THEN 00710 CALL rescale_vel_region(part,md_env,simpar,error=error) 00711 ELSE 00712 CALL rescale_vel(part,simpar,ekin,vcom=vcom) 00713 END IF 00714 00715 ! After rescaling, the core and shell velocities must also adapt. 00716 DO i = 1, natoms 00717 shell_index = part(i)%shell_index 00718 IF(shell_present .AND. shell_index/=0) THEN 00719 atomic_kind => part(i)%atomic_kind 00720 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, shell=shell) 00721 fac_masss = shell%mass_shell/mass 00722 fac_massc = shell%mass_core/mass 00723 vs = shell_part(shell_index)%v 00724 vc = core_part(shell_index)%v 00725 00726 shell_part(shell_index)%v(1) = part(i)%v(1) + fac_massc*(vs(1)-vc(1)) 00727 shell_part(shell_index)%v(2) = part(i)%v(2) + fac_massc*(vs(2)-vc(2)) 00728 shell_part(shell_index)%v(3) = part(i)%v(3) + fac_massc*(vs(3)-vc(3)) 00729 core_part(shell_index)%v(1) = part(i)%v(1) + fac_masss*(vc(1)-vs(1)) 00730 core_part(shell_index)%v(2) = part(i)%v(2) + fac_masss*(vc(2)-vs(2)) 00731 core_part(shell_index)%v(3) = part(i)%v(3) + fac_masss*(vc(3)-vs(3)) 00732 END IF 00733 END DO 00734 END IF 00735 ELSE 00736 ! Initializing velocities deterministically on all processors, if not given in input 00737 DO i = 1, natoms 00738 atomic_kind => part(i)%atomic_kind 00739 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 00740 part(i)%v(1) = 0.0_dp 00741 part(i)%v(2) = 0.0_dp 00742 part(i)%v(3) = 0.0_dp 00743 IF (mass.NE.0.0) THEN 00744 SELECT CASE(is_fixed(i)) 00745 CASE (use_perd_x) 00746 part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) 00747 part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) 00748 CASE (use_perd_y) 00749 part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) 00750 part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) 00751 CASE (use_perd_z) 00752 part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) 00753 part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) 00754 CASE (use_perd_xy) 00755 part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) 00756 CASE (use_perd_xz) 00757 part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) 00758 CASE (use_perd_yz) 00759 part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) 00760 CASE (use_perd_none) 00761 part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) 00762 part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) 00763 part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass) 00764 END SELECT 00765 END IF 00766 END DO 00767 00768 ! Subtract the vcom 00769 CALL compute_vcom(part,is_fixed,vcom) 00770 CALL subtract_vcom(part,is_fixed,vcom) 00771 ! If requested and the system is not periodic, subtract the angular velocity 00772 CALL force_env_get(force_env, cell=cell, error=error) 00773 IF (SUM(cell%perd(1:3)) == 0 .AND. simpar%angvel_zero) THEN 00774 CALL compute_rcom(part,is_fixed,rcom) 00775 CALL compute_vang(part,is_fixed,rcom,vang,error) 00776 CALL subtract_vang(part,is_fixed,rcom,vang) 00777 END IF 00778 ! Rescale the velocities 00779 IF(simpar%do_thermal_region) THEN 00780 CALL rescale_vel_region(part,md_env,simpar,error=error) 00781 ELSE 00782 ekin = compute_ekin(part) 00783 CALL rescale_vel(part,simpar,ekin) 00784 END IF 00785 00786 ! Initialize the core and the shell velocity. Atom velocities are just 00787 ! copied so that the initial relative core-shell velocity is zero. 00788 IF (shell_present) THEN 00789 CALL optimize_shell_core(force_env,part,shell_part, core_part, globenv, error=error) 00790 ENDIF 00791 END IF 00792 00793 IF (iw>0) THEN 00794 ! Recompute vcom, ecom and ekin for IO 00795 CALL compute_vcom(part,is_fixed,vcom,ecom) 00796 ekin = compute_ekin(part) - ecom 00797 IF (simpar%nfree == 0) THEN 00798 CPPostcondition(ekin==0.0_dp,cp_failure_level,routineP,error,failure) 00799 temp = 0.0_dp 00800 ELSE 00801 temp = 2.0_dp * ekin / REAL ( simpar%nfree,KIND=dp) 00802 END IF 00803 tmp_r1 = cp_unit_from_cp2k(temp,"K",error=error) 00804 WRITE (iw, '( A, T61, F18.2, A2 )' ) ' Initial Temperature ', tmp_r1, " K" 00805 WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )' ) ' COM velocity:', vcom ( 1 ), vcom ( 2 ), vcom ( 3 ) 00806 00807 ! compute and log rcom and vang if not periodic 00808 CALL force_env_get(force_env, cell=cell, error=error) 00809 IF (SUM(cell%perd(1:3)) == 0) THEN 00810 CALL compute_rcom(part,is_fixed,rcom) 00811 CALL compute_vang(part,is_fixed,rcom,vang,error) 00812 WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )' ) ' COM position:', rcom ( 1 ), rcom ( 2 ), rcom ( 3 ) 00813 WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )' ) ' Angular velocity:', vang ( 1 ), vang ( 2 ), vang ( 3 ) 00814 END IF 00815 WRITE ( iw, '( 1X,79("*"),/)' ) 00816 END IF 00817 00818 DEALLOCATE (is_fixed,STAT=stat) 00819 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00820 CALL cp_print_key_finished_output(iw,logger,print_section,"PROGRAM_RUN_INFO", error=error) 00821 CALL timestop(handle) 00822 00823 END SUBROUTINE initialize_velocities 00824 00825 ! ***************************************************************************** 00830 SUBROUTINE reset_vcom(subsys, md_ener, vsubtract, error) 00831 TYPE(cp_subsys_type), POINTER :: subsys 00832 TYPE(md_ener_type), POINTER :: md_ener 00833 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: vsubtract 00834 TYPE(cp_error_type), INTENT(inout) :: error 00835 00836 CHARACTER(LEN=*), PARAMETER :: routineN = 'reset_vcom', 00837 routineP = moduleN//':'//routineN 00838 00839 INTEGER :: atom, handle, iatom, ikind, 00840 natom, shell_index 00841 INTEGER, DIMENSION(:), POINTER :: atom_list 00842 LOGICAL :: is_shell 00843 REAL(KIND=dp) :: ekin_old, imass_c, imass_s, 00844 mass, v2 00845 REAL(KIND=dp), DIMENSION(3) :: tmp, v 00846 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 00847 TYPE(atomic_kind_type), POINTER :: atomic_kind 00848 TYPE(particle_list_type), POINTER :: core_particles, particles, 00849 shell_particles 00850 TYPE(shell_kind_type), POINTER :: shell 00851 00852 NULLIFY(particles, atomic_kind, atomic_kinds, atom_list, shell) 00853 CALL timeset(routineN,handle) 00854 00855 CALL cp_subsys_get(subsys,& 00856 atomic_kinds=atomic_kinds,& 00857 particles=particles,& 00858 shell_particles=shell_particles,& 00859 core_particles=core_particles,& 00860 error=error) 00861 00862 ekin_old = md_ener%ekin 00863 ! Possibly subtract a quantity from all velocities 00864 DO ikind=1,atomic_kinds%n_els 00865 atomic_kind => atomic_kinds%els(ikind) 00866 CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list,& 00867 natom=natom, mass=mass, shell_active=is_shell, shell=shell) 00868 IF (is_shell) THEN 00869 tmp = 0.5_dp*vsubtract*mass 00870 imass_s = 1.0_dp/shell%mass_shell 00871 imass_c = 1.0_dp/shell%mass_core 00872 DO iatom = 1, natom 00873 atom = atom_list(iatom) 00874 shell_index = particles%els(atom)%shell_index 00875 shell_particles%els(shell_index)%v = shell_particles%els(shell_index)%v - tmp*imass_s 00876 core_particles%els(shell_index)%v = core_particles%els(shell_index)%v - tmp*imass_c 00877 particles%els(atom)%v = particles%els(atom)%v - vsubtract 00878 END DO 00879 ELSE 00880 DO iatom = 1, natom 00881 atom = atom_list(iatom) 00882 particles%els(atom)%v = particles%els(atom)%v - vsubtract 00883 END DO 00884 END IF 00885 END DO 00886 ! Compute Kinetic Energy and COM Velocity 00887 md_ener%vcom = 0.0_dp 00888 md_ener%total_mass = 0.0_dp 00889 md_ener%ekin = 0.0_dp 00890 DO ikind=1,atomic_kinds%n_els 00891 atomic_kind => atomic_kinds%els(ikind) 00892 CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass, natom=natom) 00893 v2 = 0.0_dp 00894 v = 0.0_dp 00895 DO iatom = 1, natom 00896 atom = atom_list(iatom) 00897 v2 = v2 + SUM(particles%els(atom)%v**2) 00898 v(1) = v(1) + particles%els(atom)%v(1) 00899 v(2) = v(2) + particles%els(atom)%v(2) 00900 v(3) = v(3) + particles%els(atom)%v(3) 00901 END DO 00902 md_ener%ekin = md_ener%ekin + 0.5_dp*mass*v2 00903 md_ener%vcom(1) = md_ener%vcom(1) + mass*v(1) 00904 md_ener%vcom(2) = md_ener%vcom(2) + mass*v(2) 00905 md_ener%vcom(3) = md_ener%vcom(3) + mass*v(3) 00906 md_ener%total_mass = md_ener%total_mass + REAL(natom,KIND=dp)*mass 00907 END DO 00908 md_ener%vcom = md_ener%vcom / md_ener%total_mass 00909 md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin 00910 IF (md_ener%nfree /=0) THEN 00911 md_ener%temp_part = 2.0_dp*md_ener%ekin/REAL(md_ener%nfree,KIND=dp)*kelvin 00912 END IF 00913 CALL timestop(handle) 00914 00915 END SUBROUTINE reset_vcom 00916 00917 ! ***************************************************************************** 00922 SUBROUTINE scale_velocity(subsys, md_ener, temp_expected, temp_tol, iw, error) 00923 TYPE(cp_subsys_type), POINTER :: subsys 00924 TYPE(md_ener_type), POINTER :: md_ener 00925 REAL(KIND=dp), INTENT(IN) :: temp_expected, temp_tol 00926 INTEGER, INTENT(IN) :: iw 00927 TYPE(cp_error_type), INTENT(inout) :: error 00928 00929 CHARACTER(LEN=*), PARAMETER :: routineN = 'scale_velocity', 00930 routineP = moduleN//':'//routineN 00931 00932 REAL(KIND=dp) :: ekin_old, scale, temp_old 00933 00934 IF (ABS(temp_expected - md_ener%temp_part/kelvin) > temp_tol) THEN 00935 scale = 0.0_dp 00936 IF (md_ener%temp_part>0.0_dp) scale = SQRT((temp_expected/md_ener%temp_part)*kelvin) 00937 ekin_old = md_ener%ekin 00938 temp_old = md_ener%temp_part 00939 md_ener%ekin = 0.0_dp 00940 md_ener%temp_part = 0.0_dp 00941 md_ener%vcom = 0.0_dp 00942 md_ener%total_mass = 0.0_dp 00943 00944 CALL scale_velocity_low(subsys,scale,ireg=0,ekin=md_ener%ekin,vcom=md_ener%vcom,error=error) 00945 IF(md_ener%nfree /=0) THEN 00946 md_ener%temp_part = 2.0_dp*md_ener%ekin/REAL(md_ener%nfree,KIND=dp)*kelvin 00947 END IF 00948 md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin 00949 00950 IF (iw>0) THEN 00951 WRITE (UNIT=iw,FMT="(/,T2,A,F10.2,A,F10.2,A)")"Temperature scaled to requested temperature:",& 00952 temp_old," K ->",md_ener%temp_part," K" 00953 END IF 00954 END IF 00955 END SUBROUTINE scale_velocity 00956 00957 ! ***************************************************************************** 00961 SUBROUTINE scale_velocity_region(md_env, subsys, md_ener, simpar, iw, error) 00962 00963 TYPE(md_environment_type), POINTER :: md_env 00964 TYPE(cp_subsys_type), POINTER :: subsys 00965 TYPE(md_ener_type), POINTER :: md_ener 00966 TYPE(simpar_type), POINTER :: simpar 00967 INTEGER, INTENT(IN) :: iw 00968 TYPE(cp_error_type), INTENT(inout) :: error 00969 00970 CHARACTER(LEN=*), PARAMETER :: routineN = 'scale_velocity_region', 00971 routineP = moduleN//':'//routineN 00972 00973 INTEGER :: ireg, istat, nfree, 00974 nfree_done, nregions 00975 REAL(KIND=dp) :: ekin, ekin_old, 00976 ekin_total_new, fscale, 00977 vcom(3), vcom_total(3) 00978 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: temp_new, temp_old 00979 TYPE(particle_list_type), POINTER :: particles 00980 TYPE(particle_type), DIMENSION(:), 00981 POINTER :: part 00982 TYPE(thermal_region_type), POINTER :: t_region 00983 TYPE(thermal_regions_type), POINTER :: thermal_regions 00984 00985 NULLIFY( particles, part, thermal_regions, t_region) 00986 CALL cp_subsys_get(subsys, particles=particles, error=error) 00987 part => particles%els 00988 CALL get_md_env(md_env, thermal_regions=thermal_regions, error=error) 00989 00990 nregions = thermal_regions%nregions 00991 nfree_done = 0 00992 ekin_total_new = 0.0_dp 00993 ekin_old = md_ener%ekin 00994 vcom_total = 0.0_dp 00995 ALLOCATE(temp_new(0:nregions),temp_old(0:nregions), STAT=istat) 00996 temp_new = 0.0_dp 00997 temp_old = 0.0_dp 00998 !loop regions 00999 DO ireg = 1,nregions 01000 NULLIFY(t_region) 01001 t_region => thermal_regions%thermal_region(ireg) 01002 nfree = 3.0_dp*t_region%npart 01003 ekin = compute_ekin(part,ireg) 01004 IF(nfree > 0) t_region%temperature = 2.0_dp*ekin/REAL(nfree,KIND=dp)*kelvin 01005 temp_old(ireg) = t_region%temperature 01006 IF(t_region%temp_tol > 0.0_dp .AND. & 01007 ABS(t_region%temp_expected - t_region%temperature/kelvin) > t_region%temp_tol) THEN 01008 fscale = SQRT((t_region%temp_expected/t_region%temperature)*kelvin) 01009 CALL scale_velocity_low(subsys,fscale,ireg,ekin,vcom,error) 01010 t_region%temperature = 2.0_dp*ekin/REAL(nfree,KIND=dp)*kelvin 01011 temp_new(ireg) = t_region%temperature 01012 END IF 01013 nfree_done = nfree_done + nfree 01014 ekin_total_new = ekin_total_new + ekin 01015 END DO 01016 nfree = simpar%nfree - nfree_done 01017 ekin = compute_ekin(part,ireg=0) 01018 IF(nfree>0) thermal_regions%temp_reg0 = 2.0_dp*ekin/REAL(nfree,KIND=dp)*kelvin 01019 temp_old(0) = thermal_regions%temp_reg0 01020 IF(simpar%temp_tol > 0.0_dp .AND. nfree>0) THEN 01021 IF (ABS(simpar%temp_ext - thermal_regions%temp_reg0/kelvin) > simpar%temp_tol) THEN 01022 fscale = SQRT((simpar%temp_ext/thermal_regions%temp_reg0)*kelvin) 01023 CALL scale_velocity_low(subsys,fscale,0,ekin,vcom,error) 01024 thermal_regions%temp_reg0 = 2.0_dp*ekin/REAL(nfree,KIND=dp)*kelvin 01025 temp_new(0) = thermal_regions%temp_reg0 01026 END IF 01027 END IF 01028 ekin_total_new = ekin_total_new + ekin 01029 01030 md_ener%ekin = ekin_total_new 01031 IF(md_ener%nfree /=0) THEN 01032 md_ener%temp_part = 2.0_dp*md_ener%ekin/REAL(md_ener%nfree,KIND=dp)*kelvin 01033 END IF 01034 md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin 01035 01036 DO ireg = 0,nregions 01037 IF (iw>0 .AND. temp_new(ireg)>0.0_dp) THEN 01038 WRITE (UNIT=iw,FMT="(/,T2,A,I5, A,F10.2,A,F10.2,A)")"Temperature region ", ireg, & 01039 " rescaled from:",temp_old(ireg)," K to ",temp_new(ireg)," K" 01040 END IF 01041 END DO 01042 DEALLOCATE(temp_new,temp_old,STAT=istat) 01043 01044 END SUBROUTINE scale_velocity_region 01045 01046 ! ***************************************************************************** 01050 SUBROUTINE scale_velocity_low(subsys,fscale,ireg,ekin,vcom,error) 01051 01052 TYPE(cp_subsys_type), POINTER :: subsys 01053 REAL(KIND=dp), INTENT(IN) :: fscale 01054 INTEGER, INTENT(IN) :: ireg 01055 REAL(KIND=dp), INTENT(OUT) :: ekin, vcom(3) 01056 TYPE(cp_error_type), INTENT(inout) :: error 01057 01058 CHARACTER(LEN=*), PARAMETER :: routineN = 'scale_velocity_low', 01059 routineP = moduleN//':'//routineN 01060 01061 INTEGER :: atom, iatom, ikind, my_ireg, 01062 natom, shell_index 01063 INTEGER, DIMENSION(:), POINTER :: atom_list 01064 LOGICAL :: is_shell 01065 REAL(KIND=dp) :: imass, mass, tmass, v2 01066 REAL(KIND=dp), DIMENSION(3) :: tmp, v, vc, vs 01067 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 01068 TYPE(atomic_kind_type), POINTER :: atomic_kind 01069 TYPE(particle_list_type), POINTER :: core_particles, particles, 01070 shell_particles 01071 TYPE(shell_kind_type), POINTER :: shell 01072 01073 NULLIFY(atomic_kinds,particles,shell_particles,core_particles,shell,atom_list) 01074 01075 my_ireg = ireg 01076 ekin = 0.0_dp 01077 tmass = 0.0_dp 01078 vcom = 0.0_dp 01079 01080 CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, particles=particles,& 01081 shell_particles=shell_particles, core_particles=core_particles, error=error) 01082 01083 DO ikind=1,atomic_kinds%n_els 01084 atomic_kind => atomic_kinds%els(ikind) 01085 CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass,& 01086 natom=natom, shell_active=is_shell, shell=shell) 01087 IF(is_shell) THEN 01088 imass = 1.0_dp / mass 01089 v2 = 0.0_dp 01090 v = 0.0_dp 01091 DO iatom=1,natom 01092 atom = atom_list(iatom) 01093 !check region 01094 IF( particles%els(atom)%t_region_index/=my_ireg) CYCLE 01095 01096 particles%els(atom)%v(:) = fscale*particles%els(atom)%v 01097 shell_index = particles%els(atom)%shell_index 01098 vs = shell_particles%els(shell_index)%v 01099 vc = core_particles %els(shell_index)%v 01100 tmp(1) = imass*(vs(1)-vc(1)) 01101 tmp(2) = imass*(vs(2)-vc(2)) 01102 tmp(3) = imass*(vs(3)-vc(3)) 01103 01104 shell_particles%els(shell_index)%v(1) = particles%els(atom)%v(1) + tmp(1)*shell%mass_core 01105 shell_particles%els(shell_index)%v(2) = particles%els(atom)%v(2) + tmp(2)*shell%mass_core 01106 shell_particles%els(shell_index)%v(3) = particles%els(atom)%v(3) + tmp(3)*shell%mass_core 01107 01108 core_particles%els(shell_index)%v(1) = particles%els(atom)%v(1) - tmp(1)*shell%mass_shell 01109 core_particles%els(shell_index)%v(2) = particles%els(atom)%v(2) - tmp(2)*shell%mass_shell 01110 core_particles%els(shell_index)%v(3) = particles%els(atom)%v(3) - tmp(3)*shell%mass_shell 01111 01112 ! kinetic energy and velocity of COM 01113 v2 = v2 + SUM(particles%els(atom)%v**2) 01114 v(1) = v(1) + particles%els(atom)%v(1) 01115 v(2) = v(2) + particles%els(atom)%v(2) 01116 v(3) = v(3) + particles%els(atom)%v(3) 01117 tmass = tmass + mass 01118 END DO 01119 ELSE 01120 v2 = 0.0_dp 01121 v = 0.0_dp 01122 DO iatom=1,natom 01123 atom = atom_list(iatom) 01124 !check region 01125 IF( particles%els(atom)%t_region_index/=my_ireg) CYCLE 01126 01127 particles%els(atom)%v(:) = fscale*particles%els(atom)%v 01128 ! kinetic energy and velocity of COM 01129 v2 = v2 + SUM(particles%els(atom)%v**2) 01130 v(1) = v(1) + particles%els(atom)%v(1) 01131 v(2) = v(2) + particles%els(atom)%v(2) 01132 v(3) = v(3) + particles%els(atom)%v(3) 01133 tmass = tmass + mass 01134 END DO 01135 END IF 01136 ekin = ekin + 0.5_dp*mass*v2 01137 vcom(1) = vcom(1) + mass*v(1) 01138 vcom(2) = vcom(2) + mass*v(2) 01139 vcom(3) = vcom(3) + mass*v(3) 01140 01141 END DO 01142 vcom = vcom / tmass 01143 01144 END SUBROUTINE scale_velocity_low 01145 01146 ! ***************************************************************************** 01151 SUBROUTINE scale_velocity_internal(subsys, md_ener, temp_expected, temp_tol, iw, error) 01152 TYPE(cp_subsys_type), POINTER :: subsys 01153 TYPE(md_ener_type), POINTER :: md_ener 01154 REAL(KIND=dp), INTENT(IN) :: temp_expected, temp_tol 01155 INTEGER, INTENT(IN) :: iw 01156 TYPE(cp_error_type), INTENT(inout) :: error 01157 01158 CHARACTER(LEN=*), PARAMETER :: routineN = 'scale_velocity_internal', 01159 routineP = moduleN//':'//routineN 01160 01161 INTEGER :: atom, iatom, ikind, natom, 01162 shell_index 01163 INTEGER, DIMENSION(:), POINTER :: atom_list 01164 LOGICAL :: is_shell 01165 REAL(KIND=dp) :: ekin_shell_old, fac_mass, 01166 mass, scale, temp_shell_old, 01167 v2 01168 REAL(KIND=dp), DIMENSION(3) :: tmp, v, vc, vs 01169 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 01170 TYPE(atomic_kind_type), POINTER :: atomic_kind 01171 TYPE(particle_list_type), POINTER :: core_particles, particles, 01172 shell_particles 01173 TYPE(shell_kind_type), POINTER :: shell 01174 01175 NULLIFY(atom_list,atomic_kinds,atomic_kind,core_particles,particles,shell_particles,shell) 01176 IF (ABS(temp_expected - md_ener%temp_shell/kelvin) > temp_tol) THEN 01177 scale = 0.0_dp 01178 IF (md_ener%temp_shell>EPSILON(0.0_dp)) scale = SQRT((temp_expected/md_ener%temp_shell)*kelvin) 01179 ekin_shell_old = md_ener%ekin_shell 01180 temp_shell_old = md_ener%temp_shell 01181 md_ener%ekin_shell = 0.0_dp 01182 md_ener%temp_shell = 0.0_dp 01183 01184 CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, particles=particles, shell_particles=shell_particles,& 01185 core_particles=core_particles, error=error) 01186 01187 DO ikind=1,atomic_kinds%n_els 01188 atomic_kind => atomic_kinds%els(ikind) 01189 CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass, natom=natom,& 01190 shell_active=is_shell, shell=shell) 01191 IF(is_shell) THEN 01192 fac_mass = 1.0_dp/mass 01193 v2 = 0.0_dp 01194 DO iatom= 1, natom 01195 atom = atom_list(iatom) 01196 shell_index = particles%els(atom)%shell_index 01197 vs = shell_particles%els(shell_index)%v 01198 vc = core_particles%els(shell_index)%v 01199 v = particles%els(atom)%v 01200 tmp(1) = fac_mass*(vc(1)-vs(1)) 01201 tmp(2) = fac_mass*(vc(2)-vs(2)) 01202 tmp(3) = fac_mass*(vc(3)-vs(3)) 01203 01204 shell_particles%els(shell_index)%v(1) = v(1) - shell%mass_core*scale*tmp(1) 01205 shell_particles%els(shell_index)%v(2) = v(2) - shell%mass_core*scale*tmp(2) 01206 shell_particles%els(shell_index)%v(3) = v(3) - shell%mass_core*scale*tmp(3) 01207 01208 core_particles%els(shell_index)%v(1) = v(1) + shell%mass_shell*scale*tmp(1) 01209 core_particles%els(shell_index)%v(2) = v(2) + shell%mass_shell*scale*tmp(2) 01210 core_particles%els(shell_index)%v(3) = v(3) + shell%mass_shell*scale*tmp(3) 01211 01212 vs = shell_particles%els(shell_index)%v 01213 vc = core_particles%els(shell_index)%v 01214 tmp(1) = vc(1) - vs(1) 01215 tmp(2) = vc(2) - vs(2) 01216 tmp(3) = vc(3) - vs(3) 01217 v2 = v2 + SUM(tmp**2) 01218 END DO 01219 md_ener%ekin_shell = md_ener%ekin_shell + 0.5_dp*shell%mass_core*shell%mass_shell*fac_mass*v2 01220 END IF 01221 END DO 01222 IF(md_ener%nfree_shell>0)THEN 01223 md_ener%temp_shell = 2.0_dp*md_ener%ekin_shell/REAL(md_ener%nfree_shell,KIND=dp)*kelvin 01224 END IF 01225 md_ener%constant = md_ener%constant - ekin_shell_old + md_ener%ekin_shell 01226 01227 IF (iw>0) THEN 01228 WRITE (UNIT=iw,FMT="(/,T2,A,F10.2,A,F10.2,A)")& 01229 "Temperature shell internal motion scaled to requested temperature:",& 01230 temp_shell_old," K ->",md_ener%temp_shell," K" 01231 END IF 01232 ENDIF 01233 END SUBROUTINE scale_velocity_internal 01234 01235 ! ***************************************************************************** 01240 SUBROUTINE scale_velocity_baro(md_env, md_ener, temp_expected, temp_tol, iw, error) 01241 TYPE(md_environment_type), POINTER :: md_env 01242 TYPE(md_ener_type), POINTER :: md_ener 01243 REAL(KIND=dp), INTENT(IN) :: temp_expected, temp_tol 01244 INTEGER, INTENT(IN) :: iw 01245 TYPE(cp_error_type), INTENT(inout) :: error 01246 01247 CHARACTER(LEN=*), PARAMETER :: routineN = 'scale_velocity_baro', 01248 routineP = moduleN//':'//routineN 01249 01250 INTEGER :: i, j, nfree 01251 REAL(KIND=dp) :: ekin_old, scale, temp_old 01252 TYPE(npt_info_type), POINTER :: npt( :, : ) 01253 TYPE(simpar_type), POINTER :: simpar 01254 01255 NULLIFY( npt, simpar) 01256 CALL get_md_env ( md_env, simpar = simpar, npt = npt, error=error) 01257 IF (ABS(temp_expected - md_ener%temp_baro/kelvin) > temp_tol) THEN 01258 scale = 0.0_dp 01259 IF (md_ener%temp_baro>0.0_dp) scale = SQRT((temp_expected/md_ener%temp_baro)*kelvin) 01260 ekin_old = md_ener%baro_kin 01261 temp_old = md_ener%temp_baro 01262 md_ener%baro_kin = 0.0_dp 01263 md_ener%temp_baro = 0.0_dp 01264 IF ( simpar%ensemble==npt_i_ensemble .OR. simpar%ensemble==npe_i_ensemble) THEN 01265 npt ( 1, 1 )%v = npt ( 1, 1 )%v*scale 01266 md_ener%baro_kin = 0.5_dp * npt ( 1, 1 )%v**2 * npt ( 1, 1 )%mass 01267 ELSE IF (simpar%ensemble==npt_f_ensemble .OR. simpar%ensemble==npe_f_ensemble) THEN 01268 md_ener%baro_kin = 0.0_dp 01269 DO i = 1, 3 01270 DO j = 1, 3 01271 npt(i,j)%v = npt(i,j)%v*scale 01272 md_ener%baro_kin = md_ener%baro_kin + 0.5_dp * npt(i,j)%v**2 * npt ( i, j )%mass 01273 END DO 01274 END DO 01275 END IF 01276 01277 nfree = SIZE ( npt, 1 ) * SIZE ( npt, 2 ) 01278 md_ener%temp_baro = 2.0_dp * md_ener%baro_kin / REAL(nfree,dp)*kelvin 01279 IF (iw>0) THEN 01280 WRITE (UNIT=iw,FMT="(/,T2,A,F10.2,A,F10.2,A)")& 01281 "Temperature of barostat motion scaled to requested temperature:",& 01282 temp_old," K ->",md_ener%temp_baro," K" 01283 END IF 01284 END IF 01285 END SUBROUTINE scale_velocity_baro 01286 01287 ! ***************************************************************************** 01294 SUBROUTINE temperature_control(simpar, md_env, md_ener,force_env,logger, error) 01295 01296 TYPE(simpar_type), POINTER :: simpar 01297 TYPE(md_environment_type), POINTER :: md_env 01298 TYPE(md_ener_type), POINTER :: md_ener 01299 TYPE(force_env_type), POINTER :: force_env 01300 TYPE(cp_logger_type), POINTER :: logger 01301 TYPE(cp_error_type), INTENT(inout) :: error 01302 01303 CHARACTER(LEN=*), PARAMETER :: routineN = 'temperature_control', 01304 routineP = moduleN//':'//routineN 01305 01306 INTEGER :: handle, iw 01307 LOGICAL :: failure 01308 TYPE(cp_para_env_type), POINTER :: para_env 01309 TYPE(cp_subsys_type), POINTER :: subsys 01310 01311 CALL timeset(routineN,handle) 01312 NULLIFY(subsys, para_env) 01313 CPPrecondition(ASSOCIATED(simpar),cp_failure_level,routineP,error,failure) 01314 CPPrecondition(ASSOCIATED(md_ener),cp_failure_level,routineP,error,failure) 01315 CPPrecondition(ASSOCIATED(force_env),cp_failure_level,routineP,error,failure) 01316 CALL force_env_get(force_env,subsys=subsys,para_env=para_env,error=error) 01317 iw = cp_print_key_unit_nr(logger,force_env%root_section,"MOTION%MD%PRINT%PROGRAM_RUN_INFO",& 01318 extension=".mdLog",error=error) 01319 01320 ! Control the particle motion 01321 IF(simpar%do_thermal_region) THEN 01322 CALL scale_velocity_region(md_env, subsys, md_ener, simpar, iw, error) 01323 ELSE 01324 IF (simpar%temp_tol > 0.0_dp ) THEN 01325 CALL scale_velocity(subsys, md_ener, simpar%temp_ext, simpar%temp_tol, iw, error) 01326 END IF 01327 END IF 01328 ! Control the internal core-shell motion 01329 IF(simpar%temp_sh_tol > 0.0_dp) THEN 01330 CALL scale_velocity_internal(subsys, md_ener, simpar%temp_sh_ext, simpar%temp_sh_tol, iw, error) 01331 END IF 01332 ! Control cell motion 01333 SELECT CASE (simpar%ensemble) 01334 CASE( nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, & 01335 npt_f_ensemble, npt_i_ensemble, npe_f_ensemble, npe_i_ensemble) 01336 IF(simpar%temp_baro_tol > 0.0_dp) THEN 01337 CALL scale_velocity_baro(md_env, md_ener, simpar%temp_baro_ext, simpar%temp_baro_tol, iw, error) 01338 END IF 01339 END SELECT 01340 01341 CALL cp_print_key_finished_output(iw,logger,force_env%root_section,& 01342 "MOTION%MD%PRINT%PROGRAM_RUN_INFO", error=error) 01343 CALL timestop(handle) 01344 END SUBROUTINE temperature_control 01345 01346 ! ***************************************************************************** 01352 SUBROUTINE comvel_control(md_ener,force_env, md_section, logger, error) 01353 01354 TYPE(md_ener_type), POINTER :: md_ener 01355 TYPE(force_env_type), POINTER :: force_env 01356 TYPE(section_vals_type), POINTER :: md_section 01357 TYPE(cp_logger_type), POINTER :: logger 01358 TYPE(cp_error_type), INTENT(inout) :: error 01359 01360 CHARACTER(LEN=*), PARAMETER :: routineN = 'comvel_control', 01361 routineP = moduleN//':'//routineN 01362 01363 INTEGER :: handle, iw 01364 LOGICAL :: explicit, failure 01365 REAL(KIND=dp) :: comvel_tol, temp_old, vel_com 01366 REAL(KIND=dp), DIMENSION(3) :: vcom_old 01367 TYPE(cp_subsys_type), POINTER :: subsys 01368 01369 CALL timeset(routineN,handle) 01370 NULLIFY(subsys) 01371 CPPrecondition(ASSOCIATED(force_env),cp_failure_level,routineP,error,failure) 01372 CALL force_env_get(force_env,subsys=subsys,error=error) 01373 01374 ! Print COMVEL and COM Position 01375 iw = cp_print_key_unit_nr(logger,md_section,"PRINT%CENTER_OF_MASS",extension=".mdLog",error=error) 01376 IF (iw>0) THEN 01377 WRITE (UNIT=iw,FMT="(/,T2,A,(T58,A3,F20.10))")& 01378 "Centre of mass motion (COM):","x =",md_ener%vcom(1),"y =",md_ener%vcom(2),"z =",md_ener%vcom(3) 01379 END IF 01380 CALL cp_print_key_finished_output(iw,logger,md_section,"PRINT%CENTER_OF_MASS", error=error) 01381 01382 ! If requested rescale COMVEL 01383 CALL section_vals_val_get(md_section,"COMVEL_TOL",explicit=explicit,error=error) 01384 IF ( explicit ) THEN 01385 CALL section_vals_val_get(md_section,"COMVEL_TOL",r_val=comvel_tol,error=error) 01386 iw = cp_print_key_unit_nr(logger,md_section,"PRINT%PROGRAM_RUN_INFO",& 01387 extension=".mdLog",error=error) 01388 vel_com = SQRT(md_ener%vcom(1)**2+md_ener%vcom(2)**2+md_ener%vcom(3)**2) 01389 01390 ! Subtract the velocity of the COM, if requested 01391 IF (vel_com > comvel_tol) THEN 01392 temp_old = md_ener%temp_part/kelvin 01393 vcom_old = md_ener%vcom 01394 CALL reset_vcom( subsys, md_ener, vsubtract=vcom_old, error=error) 01395 CALL scale_velocity(subsys, md_ener, temp_old, 0.0_dp, iw, error) 01396 IF (iw>0) THEN 01397 WRITE (UNIT=iw,FMT="(T2,'MD| ',A,3F16.10,A)") & 01398 "Old VCOM = ",vcom_old(1:3)," a.u.",& 01399 "New VCOM = ",md_ener%vcom(1:3)," a.u" 01400 END IF 01401 END IF 01402 CALL cp_print_key_finished_output(iw,logger,md_section,& 01403 "PRINT%PROGRAM_RUN_INFO", error=error) 01404 END IF 01405 01406 CALL timestop(handle) 01407 END SUBROUTINE comvel_control 01408 01409 ! ***************************************************************************** 01414 SUBROUTINE angvel_control(md_ener, force_env, md_section, logger, error) 01415 01416 TYPE(md_ener_type), POINTER :: md_ener 01417 TYPE(force_env_type), POINTER :: force_env 01418 TYPE(section_vals_type), POINTER :: md_section 01419 TYPE(cp_logger_type), POINTER :: logger 01420 TYPE(cp_error_type), INTENT(inout) :: error 01421 01422 CHARACTER(LEN=*), PARAMETER :: routineN = 'angvel_control', 01423 routineP = moduleN//':'//routineN 01424 01425 INTEGER :: handle, ifixd, 01426 imolecule_kind, iw, natoms, 01427 stat 01428 INTEGER, ALLOCATABLE, DIMENSION(:) :: is_fixed 01429 LOGICAL :: explicit, failure 01430 REAL(KIND=dp) :: angvel_tol, rcom(3), 01431 temp_old, vang(3), vang_new(3) 01432 TYPE(cell_type), POINTER :: cell 01433 TYPE(cp_subsys_type), POINTER :: subsys 01434 TYPE(fixd_constraint_type), 01435 DIMENSION(:), POINTER :: fixd_list 01436 TYPE(mol_kind_new_list_type), POINTER :: molecule_kinds 01437 TYPE(molecule_kind_type), DIMENSION(:), 01438 POINTER :: molecule_kind_set 01439 TYPE(molecule_kind_type), POINTER :: molecule_kind 01440 TYPE(particle_list_type), POINTER :: particles 01441 01442 CALL timeset(routineN,handle) 01443 ! If requested rescale ANGVEL 01444 CALL section_vals_val_get(md_section,"ANGVEL_TOL",explicit=explicit,error=error) 01445 IF ( explicit ) THEN 01446 NULLIFY(subsys, cell) 01447 CPPrecondition(ASSOCIATED(force_env),cp_failure_level,routineP,error,failure) 01448 CALL force_env_get(force_env,subsys=subsys,cell=cell,error=error) 01449 01450 IF (SUM(cell%perd(1:3)) == 0) THEN 01451 CALL section_vals_val_get(md_section,"ANGVEL_TOL",r_val=angvel_tol,error=error) 01452 iw = cp_print_key_unit_nr(logger,md_section,"PRINT%PROGRAM_RUN_INFO",& 01453 extension=".mdLog",error=error) 01454 01455 CALL cp_subsys_get(subsys,molecule_kinds_new=molecule_kinds,& 01456 particles=particles, error=error) 01457 01458 natoms = SIZE(particles%els) 01459 ! Build a list of all fixed atoms (if any) 01460 ALLOCATE (is_fixed(natoms),STAT=stat) 01461 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01462 01463 is_fixed = use_perd_none 01464 molecule_kind_set => molecule_kinds%els 01465 DO imolecule_kind=1,molecule_kinds%n_els 01466 molecule_kind => molecule_kind_set(imolecule_kind) 01467 CALL get_molecule_kind(molecule_kind=molecule_kind,fixd_list=fixd_list) 01468 IF (ASSOCIATED(fixd_list)) THEN 01469 DO ifixd=1,SIZE(fixd_list) 01470 IF (.NOT.fixd_list(ifixd)%restraint%active) & 01471 is_fixed(fixd_list(ifixd)%fixd) = fixd_list(ifixd)%itype 01472 END DO 01473 END IF 01474 END DO 01475 01476 ! If requested and the system is not periodic, subtract the angular velocity 01477 CALL compute_rcom(particles%els,is_fixed,rcom) 01478 CALL compute_vang(particles%els,is_fixed,rcom,vang,error) 01479 IF (SQRT(DOT_PRODUCT(vang,vang))>angvel_tol) THEN 01480 CALL subtract_vang(particles%els,is_fixed,rcom,vang) 01481 01482 ! Rescale velocities after removal 01483 temp_old = md_ener%temp_part/kelvin 01484 CALL scale_velocity(subsys, md_ener, temp_old, 0.0_dp, iw, error) 01485 CALL compute_vang(particles%els,is_fixed,rcom,vang_new,error) 01486 IF (iw>0) THEN 01487 WRITE (UNIT=iw,FMT="(T2,'MD| ',A,3F16.10,A)") & 01488 "Old VANG = ",vang(1:3)," a.u.",& 01489 "New VANG = ",vang_new(1:3)," a.u" 01490 END IF 01491 END IF 01492 01493 DEALLOCATE (is_fixed,STAT=stat) 01494 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01495 01496 CALL cp_print_key_finished_output(iw,logger,md_section,& 01497 "PRINT%PROGRAM_RUN_INFO", error=error) 01498 END IF 01499 END IF 01500 01501 CALL timestop(handle) 01502 END SUBROUTINE angvel_control 01503 01504 ! ***************************************************************************** 01509 SUBROUTINE setup_velocities(force_env, simpar, globenv, md_env, md_section, & 01510 constraint_section, write_binary_restart_file, & 01511 error) 01512 01513 TYPE(force_env_type), POINTER :: force_env 01514 TYPE(simpar_type), POINTER :: simpar 01515 TYPE(global_environment_type), POINTER :: globenv 01516 TYPE(md_environment_type), POINTER :: md_env 01517 TYPE(section_vals_type), POINTER :: md_section, constraint_section 01518 LOGICAL, INTENT(IN) :: write_binary_restart_file 01519 TYPE(cp_error_type), INTENT(inout) :: error 01520 01521 CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_velocities', 01522 routineP = moduleN//':'//routineN 01523 01524 INTEGER :: handle, nconstraint, 01525 nconstraint_fixd 01526 LOGICAL :: apply_cns0, failure, 01527 shell_adiabatic, shell_present 01528 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 01529 TYPE(cell_type), POINTER :: cell 01530 TYPE(cp_para_env_type), POINTER :: para_env 01531 TYPE(cp_subsys_type), POINTER :: subsys 01532 TYPE(mol_kind_new_list_type), POINTER :: molecule_kinds 01533 TYPE(particle_list_type), POINTER :: core_particles, particles, 01534 shell_particles 01535 TYPE(particle_type), DIMENSION(:), 01536 POINTER :: core_particle_set, 01537 particle_set, 01538 shell_particle_set 01539 TYPE(section_vals_type), POINTER :: force_env_section, 01540 print_section, subsys_section 01541 01542 failure = .FALSE. 01543 01544 CALL timeset(routineN,handle) 01545 01546 NULLIFY (atomic_kinds,cell,para_env,subsys,molecule_kinds,core_particles,particles) 01547 NULLIFY (shell_particles,core_particle_set,particle_set,shell_particle_set) 01548 NULLIFY (force_env_section,print_section,subsys_section) 01549 01550 print_section => section_vals_get_subs_vals(md_section,"PRINT",error=error) 01551 apply_cns0 = .FALSE. 01552 IF (simpar%constraint) THEN 01553 CALL section_vals_val_get(constraint_section,"CONSTRAINT_INIT",l_val=apply_cns0,error=error) 01554 END IF 01555 ! Always initialize velocities and possibly restart them 01556 CALL force_env_get(force_env, subsys=subsys, cell=cell, para_env=para_env,& 01557 force_env_section=force_env_section, error=error ) 01558 subsys_section => section_vals_get_subs_vals(force_env_section,"SUBSYS",error=error) 01559 01560 CALL cp_subsys_get(subsys,& 01561 atomic_kinds=atomic_kinds,& 01562 core_particles=core_particles,& 01563 molecule_kinds_new=molecule_kinds,& 01564 particles=particles,& 01565 shell_particles=shell_particles,& 01566 error=error) 01567 01568 CALL get_atomic_kind_set(atomic_kind_set=atomic_kinds%els,& 01569 shell_present=shell_present,& 01570 shell_adiabatic=shell_adiabatic) 01571 01572 NULLIFY (core_particle_set) 01573 NULLIFY (particle_set) 01574 NULLIFY (shell_particle_set) 01575 particle_set => particles%els 01576 01577 IF (shell_present.AND.shell_adiabatic) THEN 01578 ! Constraints are not yet implemented for core-shell models generally 01579 CALL get_molecule_kind_set(molecule_kind_set=molecule_kinds%els,& 01580 nconstraint=nconstraint,& 01581 nconstraint_fixd=nconstraint_fixd) 01582 CALL cp_assert((nconstraint - nconstraint_fixd == 0),& 01583 cp_fatal_level,cp_assertion_failed,routineP,& 01584 "Only the fixed atom constraint is implemented for core-shell models",& 01585 only_ionode=.TRUE.) 01586 !MK CPPostcondition(.NOT.simpar%constraint,cp_failure_level,routineP,error,failure) 01587 CPPostcondition(ASSOCIATED(shell_particles),cp_failure_level,routineP,error,failure) 01588 CPPostcondition(ASSOCIATED(core_particles),cp_failure_level,routineP,error,failure) 01589 shell_particle_set => shell_particles%els 01590 core_particle_set => core_particles%els 01591 END IF 01592 01593 CALL initialize_velocities(simpar,particle_set, molecule_kinds=molecule_kinds,& 01594 force_env=force_env, globenv=globenv, md_env=md_env,& 01595 label="Velocities initialization",print_section=print_section, & 01596 subsys_section=subsys_section, shell_present=(shell_present.AND.shell_adiabatic), & 01597 shell_part=shell_particle_set, core_part=core_particle_set, force_rescaling=.FALSE., & 01598 para_env=para_env,error=error) 01599 01600 ! Apply constraints if required and rescale velocities.. 01601 IF (simpar%ensemble /= reftraj_ensemble) THEN 01602 IF (apply_cns0) THEN 01603 CALL force_env_calc_energy_force ( force_env, calc_force=.TRUE.,error=error) 01604 CALL force_env_shake(force_env,shake_tol=simpar%shake_tol,& 01605 log_unit=simpar%info_constraint,lagrange_mult=simpar%lagrange_multipliers,& 01606 dump_lm=simpar%dump_lm,compold=.TRUE.,error=error) 01607 CALL force_env_rattle(force_env,shake_tol=simpar%shake_tol,& 01608 log_unit=simpar%info_constraint,lagrange_mult=simpar%lagrange_multipliers,& 01609 dump_lm=simpar%dump_lm,reset=.TRUE.,error=error) 01610 IF (simpar%do_respa)THEN 01611 CALL force_env_calc_energy_force (force_env%sub_force_env(1)%force_env,& 01612 calc_force=.TRUE.,error=error) 01613 CALL force_env_shake(force_env%sub_force_env(1)%force_env,& 01614 shake_tol=simpar%shake_tol,log_unit=simpar%info_constraint,& 01615 lagrange_mult=simpar%lagrange_multipliers,dump_lm=simpar%dump_lm,compold=.TRUE.,error=error) 01616 CALL force_env_rattle(force_env%sub_force_env(1)%force_env,& 01617 shake_tol=simpar%shake_tol,log_unit=simpar%info_constraint,& 01618 lagrange_mult=simpar%lagrange_multipliers,dump_lm=simpar%dump_lm,reset=.TRUE.,error=error) 01619 END IF 01620 ! Reinitialize velocities rescaling properly after rattle 01621 subsys_section => section_vals_get_subs_vals(force_env_section,"SUBSYS",error=error) 01622 CALL update_subsys(subsys_section,force_env,.FALSE.,write_binary_restart_file,error) 01623 01624 CALL initialize_velocities(simpar,particle_set, molecule_kinds=molecule_kinds,& 01625 force_env=force_env, globenv=globenv, md_env=md_env,& 01626 label="Re-Initializing velocities after applying constraints",print_section=print_section, & 01627 subsys_section=subsys_section, shell_present=(shell_present.AND.shell_adiabatic), & 01628 shell_part=shell_particle_set, core_part=core_particle_set, force_rescaling=.TRUE., & 01629 para_env=para_env,error=error) 01630 END IF 01631 END IF 01632 01633 ! Perform setup for a cascade run 01634 CALL initialize_cascade(simpar,particle_set,molecule_kinds,md_section,& 01635 error=error) 01636 01637 CALL timestop(handle) 01638 01639 END SUBROUTINE setup_velocities 01640 01641 ! ***************************************************************************** 01647 SUBROUTINE initialize_cascade(simpar,particle_set,molecule_kinds,md_section,& 01648 error) 01649 01650 TYPE(simpar_type), POINTER :: simpar 01651 TYPE(particle_type), DIMENSION(:), 01652 POINTER :: particle_set 01653 TYPE(mol_kind_new_list_type), POINTER :: molecule_kinds 01654 TYPE(section_vals_type), POINTER :: md_section 01655 TYPE(cp_error_type), INTENT(INOUT) :: error 01656 01657 CHARACTER(LEN=*), PARAMETER :: routineN = 'initialize_cascade', 01658 routineP = moduleN//':'//routineN 01659 01660 CHARACTER(len=2*default_string_length) :: line 01661 INTEGER :: handle, iatom, ifixd, 01662 imolecule_kind, iparticle, 01663 iw, natom, nparticle, stat 01664 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_index, is_fixed 01665 LOGICAL :: failure, init_cascade, is_ok, 01666 no_read_error 01667 REAL(KIND=dp) :: ecom, ekin, energy, norm, 01668 temp, temperature, v 01669 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: matom, weight 01670 REAL(KIND=dp), ALLOCATABLE, 01671 DIMENSION(:, :) :: vatom 01672 REAL(KIND=dp), DIMENSION(3) :: vcom 01673 TYPE(atomic_kind_type), POINTER :: atomic_kind 01674 TYPE(cp_logger_type), POINTER :: logger 01675 TYPE(cp_sll_val_type), POINTER :: atom_list 01676 TYPE(fixd_constraint_type), 01677 DIMENSION(:), POINTER :: fixd_list 01678 TYPE(molecule_kind_type), DIMENSION(:), 01679 POINTER :: molecule_kind_set 01680 TYPE(molecule_kind_type), POINTER :: molecule_kind 01681 TYPE(section_vals_type), POINTER :: atom_list_section, 01682 cascade_section, print_section 01683 TYPE(val_type), POINTER :: val 01684 01685 CALL timeset(routineN,handle) 01686 01687 failure = .FALSE. 01688 01689 NULLIFY (atom_list) 01690 NULLIFY (atom_list_section) 01691 NULLIFY (atomic_kind) 01692 NULLIFY (cascade_section) 01693 NULLIFY (fixd_list) 01694 NULLIFY (molecule_kind) 01695 NULLIFY (molecule_kind_set) 01696 NULLIFY (logger) 01697 NULLIFY (val) 01698 01699 logger => cp_error_get_logger(error) 01700 print_section => section_vals_get_subs_vals(md_section,"PRINT",error=error) 01701 iw = cp_print_key_unit_nr(logger,print_section,"PROGRAM_RUN_INFO",extension=".log",error=error) 01702 01703 cascade_section => section_vals_get_subs_vals(md_section,"CASCADE",error=error) 01704 CALL section_vals_val_get(cascade_section,"_SECTION_PARAMETERS_",l_val=init_cascade,error=error) 01705 01706 nparticle = SIZE(particle_set) 01707 01708 IF (init_cascade) THEN 01709 01710 CALL section_vals_val_get(cascade_section,"ENERGY",r_val=energy,error=error) 01711 CALL cp_assert((energy >= 0.0_dp),cp_fatal_level,cp_assertion_failed,routineP,& 01712 "Error occurred reading &CASCADE section: Negative energy found",& 01713 only_ionode=.TRUE.) 01714 01715 IF (iw > 0) THEN 01716 ekin = cp_unit_from_cp2k(energy,"keV",error=error) 01717 WRITE (UNIT=iw,FMT="(T2,A,T61,F20.6)")& 01718 "CASCADE| Energy [keV]",ekin 01719 WRITE (UNIT=iw,FMT="(T2,A)")& 01720 "CASCADE|" 01721 END IF 01722 01723 ! Read the atomic velocities given in the input file 01724 atom_list_section => section_vals_get_subs_vals(cascade_section,"ATOM_LIST",error=error) 01725 CALL section_vals_val_get(atom_list_section,"_DEFAULT_KEYWORD_",n_rep_val=natom,error=error) 01726 CALL section_vals_list_get(atom_list_section,"_DEFAULT_KEYWORD_",list=atom_list,error=error) 01727 CALL cp_assert((natom > 0),cp_fatal_level,cp_assertion_failed,routineP,& 01728 "Error occurred reading &CASCADE section: No atom list found",& 01729 only_ionode=.TRUE.) 01730 01731 IF (iw > 0) THEN 01732 WRITE (UNIT=iw,FMT="(T2,A,T11,A,3(11X,A),9X,A)")& 01733 "CASCADE| ","Atom index","v(x)","v(y)","v(z)","weight" 01734 END IF 01735 01736 ALLOCATE (atom_index(natom),STAT=stat) 01737 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01738 ALLOCATE (matom(natom),STAT=stat) 01739 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01740 ALLOCATE (vatom(3,natom),STAT=stat) 01741 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01742 ALLOCATE (weight(natom),STAT=stat) 01743 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01744 01745 DO iatom=1,natom 01746 is_ok = cp_sll_val_next(atom_list,val,error=error) 01747 CALL val_get(val,c_val=line,error=error) 01748 ! Read atomic index, velocity vector, and weight 01749 no_read_error = .FALSE. 01750 READ (UNIT=line,FMT=*,ERR=999) atom_index(iatom),vatom(1:3,iatom),weight(iatom) 01751 no_read_error = .TRUE. 01752 999 CALL cp_assert(no_read_error,cp_fatal_level,cp_assertion_failed,routineP,& 01753 "Error occurred reading &CASCADE section. Last line read <"//& 01754 TRIM(line)//">",only_ionode=.TRUE.) 01755 CALL cp_assert((atom_index(iatom) > 0).AND.((atom_index(iatom) <= nparticle)),& 01756 cp_fatal_level,cp_assertion_failed,routineP,& 01757 "Error occurred reading &CASCADE section: Invalid atom index found",& 01758 only_ionode=.TRUE.) 01759 CALL cp_assert((weight(iatom) >= 0.0_dp),cp_fatal_level,cp_assertion_failed,routineP,& 01760 "Error occurred reading &CASCADE section: Negative weight found",& 01761 only_ionode=.TRUE.) 01762 IF (iw > 0) THEN 01763 WRITE (UNIT=iw,FMT="(T2,A,I10,4(1X,F14.6))")& 01764 "CASCADE| ",atom_index(iatom),vatom(1:3,iatom),weight(iatom) 01765 END IF 01766 END DO 01767 01768 ! Normalise velocities and weights 01769 norm = 0.0_dp 01770 DO iatom=1,natom 01771 iparticle = atom_index(iatom) 01772 CALL cp_assert((particle_set(iparticle)%shell_index == 0),& 01773 cp_warning_level,cp_assertion_failed,routineP,& 01774 "Warning: The primary knock-on atom is a core-shell atom",& 01775 only_ionode=.TRUE.) 01776 atomic_kind => particle_set(iparticle)%atomic_kind 01777 CALL get_atomic_kind(atomic_kind=atomic_kind,mass=matom(iatom)) 01778 norm = norm + matom(iatom)*weight(iatom) 01779 END DO 01780 weight(:) = matom(:)*weight(:)*energy/norm 01781 DO iatom=1,natom 01782 norm = SQRT(DOT_PRODUCT(vatom(1:3,iatom),vatom(1:3,iatom))) 01783 vatom(1:3,iatom) = vatom(1:3,iatom)/norm 01784 END DO 01785 01786 IF (iw > 0) THEN 01787 WRITE (UNIT=iw,FMT="(T2,A)")& 01788 "CASCADE|",& 01789 "CASCADE| Normalised velocities and additional kinetic energy [keV]",& 01790 "CASCADE|" 01791 WRITE (UNIT=iw,FMT="(T2,A,T11,A,3(11X,A),9X,A)")& 01792 "CASCADE| ","Atom index","v(x)","v(y)","v(z)","E(kin)" 01793 DO iatom=1,natom 01794 ekin = cp_unit_from_cp2k(weight(iatom),"keV",error=error) 01795 WRITE (UNIT=iw,FMT="(T2,A,I10,4(1X,F14.6))")& 01796 "CASCADE| ",atom_index(iatom),vatom(1:3,iatom),ekin 01797 END DO 01798 END IF 01799 01800 ! Apply velocity modifications 01801 DO iatom=1,natom 01802 iparticle = atom_index(iatom) 01803 particle_set(iparticle)%v(:) = particle_set(iparticle)%v(:) +& 01804 SQRT(2.0_dp*weight(iatom)/matom(iatom))*vatom(1:3,iatom) 01805 END DO 01806 01807 DEALLOCATE (atom_index,STAT=stat) 01808 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01809 DEALLOCATE (matom,STAT=stat) 01810 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01811 DEALLOCATE (vatom,STAT=stat) 01812 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01813 DEALLOCATE (weight,STAT=stat) 01814 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01815 01816 IF (iw > 0) THEN 01817 ! Build a list of all fixed atoms (if any) 01818 ALLOCATE (is_fixed(nparticle),STAT=stat) 01819 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01820 is_fixed = use_perd_none 01821 molecule_kind_set => molecule_kinds%els 01822 DO imolecule_kind=1,molecule_kinds%n_els 01823 molecule_kind => molecule_kind_set(imolecule_kind) 01824 CALL get_molecule_kind(molecule_kind=molecule_kind,fixd_list=fixd_list) 01825 IF (ASSOCIATED(fixd_list)) THEN 01826 DO ifixd=1,SIZE(fixd_list) 01827 IF (.NOT.fixd_list(ifixd)%restraint%active) is_fixed(fixd_list(ifixd)%fixd) = fixd_list(ifixd)%itype 01828 END DO 01829 END IF 01830 END DO 01831 ! Compute vcom, ecom and ekin for printout 01832 CALL compute_vcom(particle_set,is_fixed,vcom,ecom) 01833 ekin = compute_ekin(particle_set) - ecom 01834 IF (simpar%nfree == 0) THEN 01835 CPPostcondition((ekin == 0.0_dp),cp_failure_level,routineP,error,failure) 01836 temp = 0.0_dp 01837 ELSE 01838 temp = 2.0_dp*ekin/REAL(simpar%nfree,KIND=dp) 01839 END IF 01840 temperature = cp_unit_from_cp2k(temp,"K",error=error) 01841 WRITE (UNIT=iw,FMT="(T2,A)")& 01842 "CASCADE|" 01843 WRITE (UNIT=iw,FMT="(T2,A,T61,F18.2,A2)")& 01844 "CASCADE| Temperature after cascade initialization",temperature," K" 01845 WRITE (UNIT=iw,FMT="(T2,A,T30,3(1X,ES16.8),/)")& 01846 "CASCADE| COM velocity: ",vcom(1:3) 01847 !MK ! compute and log rcom and vang if not periodic 01848 !MK CALL force_env_get(force_env,cell=cell,error=error) 01849 !MK IF (SUM(cell%perd(1:3)) == 0) THEN 01850 !MK CALL compute_rcom(particle_set,is_fixed,rcom) 01851 !MK CALL compute_vang(particle_set,is_fixed,rcom,vang,error) 01852 !MK WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )' ) ' COM position:',rcom(1:3) 01853 !MK WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )' ) ' Angular velocity:',vang(1:3) 01854 !MK END IF 01855 DEALLOCATE (is_fixed,STAT=stat) 01856 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01857 END IF 01858 01859 END IF 01860 01861 CALL cp_print_key_finished_output(iw,logger,print_section,"PROGRAM_RUN_INFO",error=error) 01862 01863 CALL timestop(handle) 01864 01865 END SUBROUTINE initialize_cascade 01866 01867 END MODULE md_vel_utils
1.7.3