|
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 ! ***************************************************************************** 00014 MODULE md_energies 00015 USE atomic_kind_list_types, ONLY: atomic_kind_list_type 00016 USE atomic_kind_types, ONLY: atomic_kind_type,& 00017 get_atomic_kind_set 00018 USE averages_types, ONLY: average_quantities_type,& 00019 compute_averages 00020 USE barostat_types, ONLY: barostat_type 00021 USE barostat_utils, ONLY: print_barostat_status 00022 USE cell_types, ONLY: cell_type,& 00023 get_cell 00024 USE cp_output_handling, ONLY: cp_p_file,& 00025 cp_print_key_finished_output,& 00026 cp_print_key_should_output,& 00027 cp_print_key_unit_nr 00028 USE cp_para_types, ONLY: cp_para_env_type 00029 USE cp_subsys_types, ONLY: cp_subsys_get,& 00030 cp_subsys_type 00031 USE cp_units, ONLY: cp_unit_from_cp2k 00032 USE f77_blas 00033 USE force_env_types, ONLY: force_env_get,& 00034 force_env_type,& 00035 use_mixed_force 00036 USE input_constants, ONLY: npe_f_ensemble,& 00037 npe_i_ensemble,& 00038 nph_uniaxial_damped_ensemble,& 00039 nph_uniaxial_ensemble,& 00040 npt_f_ensemble,& 00041 npt_i_ensemble,& 00042 reftraj_ensemble 00043 USE input_cp2k_motion, ONLY: create_md_section 00044 USE input_enumeration_types, ONLY: enum_i2c,& 00045 enumeration_type 00046 USE input_keyword_types, ONLY: keyword_get,& 00047 keyword_type 00048 USE input_section_types, ONLY: section_get_keyword,& 00049 section_release,& 00050 section_type,& 00051 section_vals_get_subs_vals,& 00052 section_vals_type,& 00053 section_vals_val_get 00054 USE kinds, ONLY: default_string_length,& 00055 dp 00056 USE machine, ONLY: m_flush 00057 USE md_conserved_quantities, ONLY: compute_conserved_quantity 00058 USE md_ener_types, ONLY: md_ener_type,& 00059 zero_md_ener 00060 USE md_environment_types, ONLY: get_md_env,& 00061 md_environment_type,& 00062 set_md_env 00063 USE motion_utils, ONLY: write_simulation_cell,& 00064 write_stress_tensor,& 00065 write_trajectory 00066 USE particle_list_types, ONLY: particle_list_type 00067 USE particle_types, ONLY: write_structure_data 00068 USE physcon, ONLY: angstrom,& 00069 femtoseconds,& 00070 kelvin 00071 USE qmmm_types, ONLY: force_mixing_label_QM_dynamics,& 00072 qmmm_env_qm_type 00073 USE reftraj_types, ONLY: reftraj_type 00074 USE simpar_types, ONLY: simpar_type 00075 USE thermal_region_types, ONLY: thermal_regions_type 00076 USE thermal_region_utils, ONLY: print_thermal_regions 00077 USE thermostat_types, ONLY: thermostats_type 00078 USE thermostat_utils, ONLY: print_thermostats_status 00079 USE timings, ONLY: timeset,& 00080 timestop 00081 USE virial_types, ONLY: virial_type 00082 #include "cp_common_uses.h" 00083 00084 IMPLICIT NONE 00085 00086 PRIVATE 00087 00088 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'md_energies' 00089 00090 00091 PUBLIC :: initialize_md_ener,& 00092 md_energy,& 00093 md_ener_reftraj,& 00094 update_nfree_qm,& 00095 md_write_output 00096 00097 CONTAINS 00098 00099 ! ***************************************************************************** 00104 SUBROUTINE initialize_md_ener(md_ener,force_env,simpar,error) 00105 00106 TYPE(md_ener_type), POINTER :: md_ener 00107 TYPE(force_env_type), POINTER :: force_env 00108 TYPE(simpar_type), POINTER :: simpar 00109 TYPE(cp_error_type), INTENT(inout) :: error 00110 00111 CHARACTER(len=*), PARAMETER :: routineN = 'initialize_md_ener', 00112 routineP = moduleN//':'//routineN 00113 00114 INTEGER :: istat, nkind 00115 LOGICAL :: failure, shell_adiabatic 00116 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 00117 TYPE(atomic_kind_type), DIMENSION(:), 00118 POINTER :: atomic_kind_set 00119 TYPE(cp_subsys_type), POINTER :: subsys 00120 TYPE(particle_list_type), POINTER :: particles, shell_particles 00121 00122 failure =.FALSE. 00123 00124 NULLIFY(subsys) 00125 NULLIFY(atomic_kinds, atomic_kind_set, particles, shell_particles) 00126 00127 CPPrecondition(ASSOCIATED(md_ener),cp_failure_level,routineP,error,failure) 00128 CPPrecondition(ASSOCIATED(force_env),cp_failure_level,routineP,error,failure) 00129 00130 CALL force_env_get(force_env, subsys=subsys, error=error) 00131 CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, particles=particles,& 00132 shell_particles=shell_particles,& 00133 error=error) 00134 atomic_kind_set => atomic_kinds%els 00135 nkind = SIZE(atomic_kind_set) 00136 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 00137 shell_adiabatic=shell_adiabatic) 00138 00139 md_ener%nfree = simpar%nfree 00140 md_ener%nfree_shell = -HUGE(0) 00141 00142 IF(shell_adiabatic) THEN 00143 md_ener%nfree_shell = 3*(shell_particles%n_els) 00144 END IF 00145 CALL update_nfree_qm(md_ener, force_env, error) 00146 00147 IF(simpar%temperature_per_kind) THEN 00148 ALLOCATE(md_ener%temp_kind(nkind), STAT=istat) 00149 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00150 ALLOCATE(md_ener%ekin_kind(nkind), STAT=istat) 00151 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00152 ALLOCATE(md_ener%nfree_kind(nkind),STAT=istat) 00153 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00154 md_ener%nfree_kind = 0 00155 00156 IF(shell_adiabatic) THEN 00157 ALLOCATE(md_ener%temp_shell_kind(nkind), STAT=istat) 00158 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00159 ALLOCATE(md_ener%ekin_shell_kind(nkind), STAT=istat) 00160 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00161 ALLOCATE(md_ener%nfree_shell_kind(nkind), STAT=istat) 00162 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00163 md_ener%nfree_shell_kind = 0 00164 END IF 00165 00166 END IF 00167 CALL zero_md_ener(md_ener, tkind=simpar%temperature_per_kind, & 00168 tshell=shell_adiabatic, error=error) 00169 md_ener%epot = 0.0_dp 00170 00171 END SUBROUTINE initialize_md_ener 00172 00173 SUBROUTINE update_nfree_qm(md_ener, force_env, error) 00174 TYPE(md_ener_type), POINTER :: md_ener 00175 TYPE(force_env_type), POINTER :: force_env 00176 TYPE(cp_error_type), INTENT(inout) :: error 00177 00178 CHARACTER(len=*), PARAMETER :: routineN = 'update_nfree_qm', 00179 routineP = moduleN//':'//routineN 00180 00181 INTEGER :: ip 00182 INTEGER, POINTER :: cur_indices(:), cur_labels(:) 00183 LOGICAL :: failure 00184 TYPE(cp_subsys_type), POINTER :: subsys 00185 TYPE(particle_list_type), POINTER :: particles 00186 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env 00187 TYPE(section_vals_type), POINTER :: force_env_section 00188 00189 failure = .FALSE. 00190 00191 md_ener%nfree_qm = -HUGE(0) 00192 00193 NULLIFY(qmmm_env,subsys,particles,force_env_section) 00194 CALL force_env_get(force_env, subsys=subsys, qmmm_env=qmmm_env, force_env_section=force_env_section, error=error) 00195 IF(ASSOCIATED(qmmm_env)) THEN 00196 IF (SIZE(force_env%sub_force_env) == 2) THEN ! doing force mixing 00197 CALL section_vals_val_get(force_env_section,"QMMM%FORCE_MIXING%RESTART_INFO%INDICES",i_vals=cur_indices,error=error) 00198 CALL section_vals_val_get(force_env_section,"QMMM%FORCE_MIXING%RESTART_INFO%LABELS",i_vals=cur_labels,error=error) 00199 md_ener%nfree_qm = 0 00200 DO ip=1, SIZE(cur_indices) 00201 IF (cur_labels(ip) >= force_mixing_label_QM_dynamics) THEN ! this is a QM atom 00202 md_ener%nfree_qm = md_ener%nfree_qm + 3 00203 END IF 00204 END DO 00205 ELSE 00206 CALL cp_subsys_get(subsys, particles=particles, error=error) 00207 ! The degrees of freedom for the quantum part of the system 00208 ! are set to 3*Number of QM atoms and to simpar%nfree in case all the MM 00209 ! system is treated at QM level (not really QM/MM, just for consistency). 00210 ! The degree of freedom will not be correct if 1-3 atoms are treated only 00211 ! MM. In this case we should take care of rotations 00212 md_ener%nfree_qm = 3*SIZE(qmmm_env%qm_atom_index) 00213 IF (md_ener%nfree_qm == 3*(particles%n_els)) md_ener%nfree_qm = md_ener%nfree 00214 ENDIF 00215 END IF 00216 00217 END SUBROUTINE update_nfree_qm 00218 00219 00220 ! ***************************************************************************** 00225 SUBROUTINE md_energy(md_env, md_ener, error) 00226 00227 TYPE(md_environment_type), POINTER :: md_env 00228 TYPE(md_ener_type), POINTER :: md_ener 00229 TYPE(cp_error_type), INTENT(inout) :: error 00230 00231 CHARACTER(len=*), PARAMETER :: routineN = 'md_energy', 00232 routineP = moduleN//':'//routineN 00233 00234 INTEGER :: natom 00235 LOGICAL :: shell_adiabatic, tkind, tshell 00236 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 00237 TYPE(atomic_kind_type), DIMENSION(:), 00238 POINTER :: atomic_kind_set 00239 TYPE(cp_subsys_type), POINTER :: subsys 00240 TYPE(force_env_type), POINTER :: force_env 00241 TYPE(particle_list_type), POINTER :: particles 00242 TYPE(simpar_type), POINTER :: simpar 00243 00244 NULLIFY(atomic_kinds, atomic_kind_set, force_env,& 00245 particles, subsys, simpar) 00246 CALL get_md_env(md_env=md_env, force_env=force_env, & 00247 simpar=simpar, error=error) 00248 00249 CALL force_env_get(force_env, & 00250 potential_energy=md_ener%epot, subsys=subsys, error=error) 00251 00252 CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, error=error) 00253 atomic_kind_set => atomic_kinds%els 00254 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 00255 shell_adiabatic=shell_adiabatic) 00256 00257 tkind = simpar%temperature_per_kind 00258 tshell = shell_adiabatic 00259 00260 CALL cp_subsys_get(subsys,particles=particles,error=error) 00261 natom=particles%n_els 00262 00263 CALL compute_conserved_quantity(md_env, md_ener, tkind=tkind,& 00264 tshell=tshell, natom=natom, error=error) 00265 00266 END SUBROUTINE md_energy 00267 00268 ! ***************************************************************************** 00275 SUBROUTINE md_ener_reftraj(md_env,md_ener,error) 00276 TYPE(md_environment_type), POINTER :: md_env 00277 TYPE(md_ener_type), POINTER :: md_ener 00278 TYPE(cp_error_type), INTENT(inout) :: error 00279 00280 CHARACTER(len=*), PARAMETER :: routineN = 'md_ener_reftraj', 00281 routineP = moduleN//':'//routineN 00282 00283 TYPE(force_env_type), POINTER :: force_env 00284 TYPE(reftraj_type), POINTER :: reftraj 00285 00286 CALL zero_md_ener(md_ener, tkind=.FALSE., tshell=.FALSE., error=error) 00287 CALL get_md_env(md_env=md_env, force_env=force_env, reftraj=reftraj, error=error) 00288 00289 IF(reftraj%info%eval_ef) THEN 00290 CALL force_env_get(force_env, potential_energy=md_ener%epot, error=error) 00291 ELSE 00292 md_ener%epot = reftraj%epot 00293 md_ener%delta_epot = (reftraj%epot - reftraj%epot0)/REAL(reftraj%natom, kind=dp)*kelvin 00294 END IF 00295 00296 END SUBROUTINE md_ener_reftraj 00297 00298 ! ***************************************************************************** 00307 SUBROUTINE md_write_output(md_env, error) 00308 00309 TYPE(md_environment_type), POINTER :: md_env 00310 TYPE(cp_error_type), INTENT(inout) :: error 00311 00312 CHARACTER(len=*), PARAMETER :: routineN = 'md_write_output', 00313 routineP = moduleN//':'//routineN 00314 00315 CHARACTER(LEN=default_string_length) :: fmd, my_act, my_pos 00316 INTEGER :: ene, ener_mix, handle, i, 00317 nat, nkind, shene, tempkind, 00318 trsl 00319 INTEGER, POINTER :: itimes 00320 LOGICAL :: failure, init, is_mixed, 00321 new_file, qmmm, 00322 shell_adiabatic, shell_present 00323 REAL(dp) :: abc( 3 ), cell_angle( 3 ), 00324 dt, econs, pv_scalar, pv_xx, 00325 pv_xx_nc 00326 REAL(KIND=dp) :: harm_shell, hugoniot 00327 REAL(KIND=dp), POINTER :: time, used_time 00328 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 00329 TYPE(atomic_kind_type), DIMENSION(:), 00330 POINTER :: atomic_kind_set 00331 TYPE(average_quantities_type), POINTER :: averages 00332 TYPE(barostat_type), POINTER :: barostat 00333 TYPE(cell_type), POINTER :: cell 00334 TYPE(cp_logger_type), POINTER :: logger 00335 TYPE(cp_para_env_type), POINTER :: para_env 00336 TYPE(cp_subsys_type), POINTER :: subsys 00337 TYPE(force_env_type), POINTER :: force_env 00338 TYPE(md_ener_type), POINTER :: md_ener 00339 TYPE(particle_list_type), POINTER :: core_particles, particles, 00340 shell_particles 00341 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env 00342 TYPE(reftraj_type), POINTER :: reftraj 00343 TYPE(section_vals_type), POINTER :: motion_section, print_key, 00344 root_section 00345 TYPE(simpar_type), POINTER :: simpar 00346 TYPE(thermal_regions_type), POINTER :: thermal_regions 00347 TYPE(thermostats_type), POINTER :: thermostats 00348 TYPE(virial_type), POINTER :: virial 00349 00350 failure = .FALSE. 00351 NULLIFY(logger) 00352 logger => cp_error_get_logger(error) 00353 CALL timeset(routineN,handle) 00354 00355 ! Zeroing 00356 hugoniot = 0.0_dp 00357 econs = 0.0_dp 00358 shell_adiabatic = .FALSE. 00359 shell_present = .FALSE. 00360 NULLIFY(motion_section, atomic_kinds, atomic_kind_set, cell, subsys, & 00361 force_env, md_ener, qmmm_env, reftraj, core_particles, particles, & 00362 shell_particles, print_key, root_section, simpar, virial, & 00363 thermostats, thermal_regions) 00364 00365 CALL get_md_env(md_env=md_env, itimes=itimes, t=time, used_time=used_time,& 00366 simpar=simpar, force_env=force_env, init=init, md_ener=md_ener,& 00367 reftraj=reftraj, thermostats=thermostats, barostat=barostat, & 00368 para_env=para_env, averages=averages, thermal_regions=thermal_regions, error=error) 00369 00370 root_section => force_env%root_section 00371 motion_section => section_vals_get_subs_vals(root_section,"MOTION",error=error) 00372 00373 CALL force_env_get(force_env, cell=cell, subsys=subsys, qmmm_env=qmmm_env,& 00374 virial=virial, error=error) 00375 00376 qmmm = (md_ener%nfree_qm>0) 00377 is_mixed = (force_env%in_use == use_mixed_force) 00378 00379 CALL cp_subsys_get(subsys,particles=particles,error=error) 00380 nat = particles%n_els 00381 dt = simpar%dt*simpar%dt_fact 00382 00383 ! Computing the scalar pressure 00384 IF ( virial%pv_availability ) THEN 00385 pv_scalar = 0._dp 00386 DO i = 1, 3 00387 pv_scalar = pv_scalar + virial%pv_total(i,i) 00388 END DO 00389 pv_scalar = pv_scalar/3._dp/cell%deth 00390 pv_scalar = cp_unit_from_cp2k(pv_scalar,"bar",error=error) 00391 pv_xx_nc = virial%pv_total(1,1)/cell%deth 00392 pv_xx = cp_unit_from_cp2k(virial%pv_total(1,1)/cell%deth,"bar",error=error) 00393 ENDIF 00394 00395 CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, error=error) 00396 atomic_kind_set => atomic_kinds%els 00397 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 00398 shell_present=shell_present, & 00399 shell_adiabatic=shell_adiabatic) 00400 00401 CALL get_cell (cell, abc=abc, alpha=cell_angle(3), beta=cell_angle(2), gamma=cell_angle(1)) 00402 00403 ! Determine POS and ACT for I/O 00404 my_pos = "APPEND" 00405 my_act = "WRITE" 00406 IF (init.AND.(itimes==0)) THEN 00407 my_pos = "REWIND" 00408 my_act = "WRITE" 00409 END IF 00410 00411 ! In case of REFTRAJ ensemble the POS is determined differently.. 00412 ! according the REFTRAJ counter 00413 IF (simpar%ensemble==reftraj_ensemble) THEN 00414 IF((reftraj%isnap==reftraj%info%first_snapshot)) THEN 00415 my_pos = "REWIND" 00416 END IF 00417 ENDIF 00418 00419 ! Performing protocol relevant to the first step of an MD run 00420 IF (init) THEN 00421 ! Computing the Hugoniot for NPH calculations 00422 IF ( simpar%ensemble == nph_uniaxial_ensemble .OR.& 00423 simpar%ensemble == nph_uniaxial_damped_ensemble ) THEN 00424 IF (simpar%e0 == 0._dp ) simpar%e0 = md_ener%epot + md_ener%ekin 00425 hugoniot = md_ener%epot + md_ener%ekin - simpar%e0 - 0.5_dp*(pv_xx_nc + simpar%p0)*& 00426 (simpar%v0-cell%deth) 00427 ENDIF 00428 00429 IF (simpar%ensemble==reftraj_ensemble) reftraj%init = init 00430 ELSE 00431 ! Performing protocol for anything beyond the first step of MD 00432 IF ( simpar%ensemble == nph_uniaxial_ensemble .OR. simpar%ensemble == nph_uniaxial_damped_ensemble) THEN 00433 hugoniot = md_ener%epot + md_ener%ekin - simpar%e0 - 0.5_dp * ( pv_xx_nc + simpar%p0 ) * & 00434 ( simpar%v0 - cell%deth ) 00435 END IF 00436 00437 IF (simpar%ensemble==reftraj_ensemble) THEN 00438 time = reftraj%time 00439 econs = md_ener%delta_epot 00440 itimes = reftraj%itimes 00441 ELSE 00442 econs = md_ener%delta_cons 00443 END IF 00444 00445 ! Compute average quantities 00446 CALL compute_averages(averages, force_env, md_ener, cell, virial, pv_scalar,& 00447 pv_xx, used_time, hugoniot, abc, cell_angle, nat, itimes, time, my_pos,& 00448 my_act, error) 00449 END IF 00450 00451 ! Print md information 00452 CALL md_write_info_low(simpar, md_ener, qmmm, virial, reftraj, cell, abc,& 00453 cell_angle, itimes, dt, time, used_time, averages, econs, pv_scalar, pv_xx,& 00454 hugoniot, nat, init, logger, motion_section, my_pos, my_act, error) 00455 00456 ! Real Ouput driven by the PRINT sections 00457 IF ((.NOT.init).OR.(itimes==0).OR.simpar%ensemble==reftraj_ensemble) THEN 00458 ! Print Energy 00459 ene = cp_print_key_unit_nr(logger,motion_section,"MD%PRINT%ENERGY",& 00460 extension=".ener",file_position=my_pos, file_action=my_act, is_new_file=new_file,& 00461 error=error) 00462 IF (ene>0) THEN 00463 IF (new_file) THEN 00464 ! Please change also the corresponding format explaination below 00465 ! keep the constant of motion the true constant of motion ! 00466 WRITE (ene,'("#",5X,A,10X,A,8X,A,10X,A,12X,A,2(8X,A))')"Step Nr.","Time[fs]","Kin.[a.u.]","Temp[K]",& 00467 "Pot.[a.u.]","Cons Qty[a.u.]","UsedTime[s]" 00468 END IF 00469 WRITE (ene,"(I10,F20.6,F20.9,F20.9,F20.9,F20.9,F20.9)")& 00470 itimes,time*femtoseconds,md_ener%ekin,md_ener%temp_part, md_ener%epot,md_ener%constant,used_time 00471 CALL m_flush(ene) 00472 END IF 00473 CALL cp_print_key_finished_output(ene,logger,motion_section,"MD%PRINT%ENERGY", error=error) 00474 00475 ! Possibly Print MIXED Energy 00476 IF (is_mixed) THEN 00477 ener_mix=cp_print_key_unit_nr(logger,motion_section,"PRINT%MIXED_ENERGIES",& 00478 extension=".ener", file_position=my_pos, file_action=my_act,& 00479 middle_name="mix", error=error) 00480 IF (ener_mix>0) THEN 00481 WRITE (ener_mix,"(I8,F12.3,F20.9,"//cp_to_string(SIZE(force_env%mixed_env%energies))//"F20.9,F20.9)")& 00482 itimes,time*femtoseconds,md_ener%epot,force_env%mixed_env%energies,md_ener%constant 00483 CALL m_flush(ener_mix) 00484 END IF 00485 CALL cp_print_key_finished_output(ener_mix,logger,motion_section,"PRINT%MIXED_ENERGIES", error=error) 00486 ENDIF 00487 00488 ! Print QMMM translation vector if requested 00489 IF (qmmm) THEN 00490 trsl = cp_print_key_unit_nr(logger,motion_section,"PRINT%TRANSLATION_VECTOR",& 00491 extension=".translation", middle_name="qmmm", error=error) 00492 IF (trsl>0) THEN 00493 WRITE(trsl,'(I10,3F15.10)')itimes,qmmm_env%transl_v 00494 END IF 00495 CALL cp_print_key_finished_output(trsl,logger,motion_section,& 00496 "PRINT%TRANSLATION_VECTOR", error=error) 00497 END IF 00498 00499 ! Write Structure data 00500 CALL write_structure_data(particles%els,cell,motion_section,error) 00501 00502 ! Print Coordinates 00503 CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot,& 00504 pos=my_pos, act=my_act, extended_xmol_title=.TRUE., error=error) 00505 00506 ! Print Velocities 00507 CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot,& 00508 "VELOCITIES", my_pos, my_act, middle_name="vel", extended_xmol_title=.TRUE., error=error) 00509 00510 ! Print Force 00511 CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot,& 00512 "FORCES", my_pos, my_act, middle_name="frc", extended_xmol_title=.TRUE., error=error) 00513 00514 ! Print Force-Mixing labels 00515 CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot,& 00516 "FORCE_MIXING_LABELS", my_pos, my_act, middle_name="fmlabels", extended_xmol_title=.TRUE., error=error) 00517 00518 ! Print Simulation Cell 00519 CALL write_simulation_cell(cell, motion_section, itimes, time*femtoseconds, my_pos, my_act, error) 00520 00521 ! Print Thermostats status 00522 CALL print_thermostats_status(thermostats, para_env, my_pos, my_act, cell, itimes, time, error) 00523 00524 ! Print Barostat status 00525 CALL print_barostat_status(barostat, simpar, my_pos, my_act, cell, itimes, time, error) 00526 00527 ! Print Stress Tensor 00528 CALL write_stress_tensor(virial, cell, motion_section, itimes, time*femtoseconds,& 00529 my_pos, my_act, error) 00530 00531 ! Temperature per Kinds 00532 IF(simpar%temperature_per_kind) THEN 00533 tempkind=cp_print_key_unit_nr(logger,motion_section,"MD%PRINT%TEMP_KIND",& 00534 extension=".temp",file_position=my_pos, file_action=my_act,error=error) 00535 IF( tempkind > 0 ) THEN 00536 nkind = SIZE(md_ener%temp_kind) 00537 fmd="(I10,F20.3,"//TRIM(ADJUSTL(cp_to_string(nkind)))//"F20.9)" 00538 fmd=TRIM(fmd) 00539 WRITE (tempkind,fmd)itimes,time*femtoseconds, md_ener%temp_kind(1:nkind) 00540 CALL m_flush(tempkind) 00541 END IF 00542 CALL cp_print_key_finished_output(tempkind,logger,motion_section,"MD%PRINT%TEMP_KIND", error=error) 00543 ELSE 00544 print_key => section_vals_get_subs_vals(motion_section,"MD%PRINT%TEMP_KIND",error=error) 00545 CALL cp_assert(.NOT.BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file),& 00546 cp_warning_level,cp_assertion_failed,routineP,& 00547 "The print_key MD%PRINT%TEMP_KIND has been activated but the "//& 00548 "calculation of the temperature per kind has not been requested. "//& 00549 "Please switch on the keyword MD%TEMP_KIND. "//& 00550 CPSourceFileRef,& 00551 only_ionode=.TRUE.) 00552 END IF 00553 !Thermal Region 00554 CALL print_thermal_regions(thermal_regions,itimes,time*femtoseconds,my_pos,my_act,error) 00555 00556 ! Core/Shell Model 00557 IF(shell_present) THEN 00558 CALL force_env_get(force_env, harmonic_shell=harm_shell, error=error) 00559 CALL cp_subsys_get(subsys, shell_particles=shell_particles, core_particles=core_particles, error=error) 00560 00561 ! Print Shell Energy 00562 shene=cp_print_key_unit_nr(logger,motion_section,"MD%PRINT%SHELL_ENERGY",& 00563 extension=".shener",file_position=my_pos, file_action=my_act, & 00564 file_form="FORMATTED", is_new_file=new_file,error=error) 00565 IF (shene>0) THEN 00566 IF(new_file) THEN 00567 WRITE (shene,'("#",3X,A,3X,A,3X,3(5X,A,5X))')"Step Nr.","Time[fs]","Kin.[a.u.]",& 00568 "Temp.[K]","Pot.[a.u.]" 00569 END IF 00570 00571 WRITE (shene,"(I8,F12.3,F20.9,F20.9,F20.9,F20.9 )")& 00572 itimes,time*femtoseconds,md_ener%ekin_shell,md_ener%temp_shell,harm_shell 00573 CALL m_flush(shene) 00574 END IF 00575 CALL cp_print_key_finished_output(shene,logger,motion_section,"MD%PRINT%SHELL_ENERGY", error=error) 00576 00577 ! Print Shell Coordinates 00578 CALL write_trajectory (force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot,& 00579 "SHELL_TRAJECTORY", my_pos, my_act, "shpos", shell_particles, extended_xmol_title=.TRUE., error=error) 00580 00581 IF(shell_adiabatic) THEN 00582 ! Print Shell Velocities 00583 CALL write_trajectory (force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot,& 00584 "SHELL_VELOCITIES", my_pos, my_act, "shvel", shell_particles, extended_xmol_title=.TRUE., error=error) 00585 00586 ! Print Shell Forces 00587 CALL write_trajectory (force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot,& 00588 "SHELL_FORCES", my_pos, my_act, "shfrc", shell_particles, extended_xmol_title=.TRUE., error=error) 00589 00590 ! Print Core Coordinates 00591 CALL write_trajectory (force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot,& 00592 "CORE_TRAJECTORY", my_pos, my_act, "copos", core_particles, extended_xmol_title=.TRUE., error=error) 00593 00594 ! Print Core Velocities 00595 CALL write_trajectory (force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot,& 00596 "CORE_VELOCITIES", my_pos, my_act, "covel", core_particles, extended_xmol_title=.TRUE., error=error) 00597 00598 ! Print Core Forces 00599 CALL write_trajectory (force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot,& 00600 "CORE_FORCES", my_pos, my_act, "cofrc", core_particles, extended_xmol_title=.TRUE., error=error) 00601 00602 ! Temperature per Kinds 00603 IF(simpar%temperature_per_kind) THEN 00604 tempkind=cp_print_key_unit_nr(logger,motion_section,"MD%PRINT%TEMP_SHELL_KIND",& 00605 extension=".shtemp", file_position=my_pos, file_action=my_act,error=error) 00606 IF( tempkind > 0 ) THEN 00607 nkind = SIZE(md_ener%temp_shell_kind) 00608 fmd="(I10,F20.3,"//TRIM(ADJUSTL(cp_to_string(nkind)))//"F20.9)" 00609 fmd=TRIM(fmd) 00610 WRITE (tempkind,fmd)itimes,time*femtoseconds, md_ener%temp_shell_kind(1:nkind) 00611 CALL m_flush(tempkind) 00612 END IF 00613 CALL cp_print_key_finished_output(tempkind, logger, motion_section,& 00614 "MD%PRINT%TEMP_SHELL_KIND", error=error) 00615 ELSE 00616 print_key => section_vals_get_subs_vals(motion_section,"MD%PRINT%TEMP_SHELL_KIND",error=error) 00617 CALL cp_assert(.NOT.BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file),& 00618 cp_warning_level,cp_assertion_failed,routineP,& 00619 "The print_key MD%PRINT%TEMP_SHELL_KIND has been activated but the "//& 00620 "calculation of the temperature per kind has not been requested. "//& 00621 "Please switch on the keyword MD%TEMP_KIND. "//& 00622 CPSourceFileRef,& 00623 only_ionode=.TRUE.) 00624 END IF 00625 END IF 00626 END IF 00627 END IF 00628 init = .FALSE. 00629 CALL set_md_env(md_env,init=init,error=error) 00630 CALL timestop(handle) 00631 END SUBROUTINE md_write_output 00632 00633 ! ***************************************************************************** 00641 SUBROUTINE md_write_info_low(simpar, md_ener, qmmm, virial, reftraj, cell,& 00642 abc, cell_angle, itimes, dt, time, used_time, averages, econs, pv_scalar, & 00643 pv_xx, hugoniot, nat, init, logger, motion_section, my_pos, my_act, error) 00644 00645 TYPE(simpar_type), POINTER :: simpar 00646 TYPE(md_ener_type), POINTER :: md_ener 00647 LOGICAL, INTENT(IN) :: qmmm 00648 TYPE(virial_type), POINTER :: virial 00649 TYPE(reftraj_type), POINTER :: reftraj 00650 TYPE(cell_type), POINTER :: cell 00651 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: abc, cell_angle 00652 INTEGER, POINTER :: itimes 00653 REAL(KIND=dp), INTENT(IN) :: dt 00654 REAL(KIND=dp), POINTER :: time, used_time 00655 TYPE(average_quantities_type), POINTER :: averages 00656 REAL(KIND=dp), INTENT(IN) :: econs, pv_scalar, pv_xx, 00657 hugoniot 00658 INTEGER, INTENT(IN) :: nat 00659 LOGICAL, INTENT(IN) :: init 00660 TYPE(cp_logger_type), POINTER :: logger 00661 TYPE(section_vals_type), POINTER :: motion_section 00662 CHARACTER(LEN=default_string_length), 00663 INTENT(IN) :: my_pos, my_act 00664 TYPE(cp_error_type), INTENT(inout) :: error 00665 00666 CHARACTER(len=*), PARAMETER :: routineN = 'md_write_info_low', 00667 routineP = moduleN//':'//routineN 00668 00669 INTEGER :: iw 00670 LOGICAL :: failure 00671 TYPE(enumeration_type), POINTER :: enum 00672 TYPE(keyword_type), POINTER :: keyword 00673 TYPE(section_type), POINTER :: section 00674 00675 failure = .FALSE. 00676 NULLIFY(enum, keyword, section) 00677 ! Print to the screen info about MD 00678 iw = cp_print_key_unit_nr(logger,motion_section,"MD%PRINT%PROGRAM_RUN_INFO",& 00679 extension=".mdLog",file_position=my_pos,file_action=my_act,error=error) 00680 CALL create_md_section(section,error=error) 00681 keyword => section_get_keyword(section,"ENSEMBLE",error=error) 00682 CALL keyword_get(keyword,enum=enum,error=error) 00683 00684 ! Performing protocol relevant to the first step of an MD run 00685 IF (iw>0) THEN 00686 IF (init) THEN 00687 ! Write initial values of quantities of interest 00688 WRITE (iw,*) 00689 WRITE (iw,'( A )') ' MD_ENERGIES| Initialization proceeding' 00690 WRITE (iw,*) 00691 WRITE (iw,'( )') 00692 WRITE (iw,'( A,A )') ' ******************************** ', & 00693 'GO CP2K GO! **********************************' 00694 WRITE (iw,'( A,A,T40,A,T60,1(1X,E20.12) )') ' INITIAL POTENTIAL ENERGY', & 00695 '[hartree]', '= ', md_ener%epot 00696 IF (simpar%ensemble/=reftraj_ensemble) THEN 00697 IF (.NOT.qmmm) THEN 00698 ! NO QM/MM info 00699 WRITE (iw,'( A,A,T40,A,T60,1(1X,E20.12) )') ' INITIAL KINETIC ENERGY', & 00700 '[hartree]', '= ', md_ener%ekin 00701 WRITE (iw,'( A,A,T40,A,T60,1(1X,F20.3) )') ' INITIAL TEMPERATURE', & 00702 '[K]', '= ', md_ener%temp_part 00703 ELSE 00704 WRITE (iw,'( A,A,T40,A,T60,1(1X,E20.12) )') ' TOTAL INITIAL KINETIC ENERGY', & 00705 '[hartree]', '= ', md_ener%ekin 00706 WRITE (iw,'( A,A,T40,A,T60,1(1X,E20.12) )') ' QM INITIAL KINETIC ENERGY', & 00707 '[hartree]', '= ', md_ener%ekin_qm 00708 WRITE (iw,'( A,A,T40,A,T60,1(1X,F20.3) )') ' TOTAL INITIAL TEMPERATURE', & 00709 '[K]', '= ', md_ener%temp_part 00710 WRITE (iw,'( A,A,T40,A,T60,1(1X,F20.3) )') ' QM INITIAL TEMPERATURE', & 00711 '[K]', '= ', md_ener%temp_qm 00712 END IF 00713 END IF 00714 IF ( simpar%ensemble == nph_uniaxial_ensemble .OR.& 00715 simpar%ensemble == nph_uniaxial_damped_ensemble.OR.& 00716 simpar%ensemble == npt_i_ensemble.OR.& 00717 simpar%ensemble == npt_f_ensemble.OR.& 00718 simpar%ensemble == npe_i_ensemble.OR.& 00719 simpar%ensemble == npe_f_ensemble) & 00720 WRITE (iw,'( A,A,T40,A,T60,1(1X,E20.12) )') ' INITIAL BAROSTAT TEMP', & 00721 '[K]', '= ', md_ener%temp_baro 00722 IF ( virial%pv_availability ) & 00723 WRITE (iw,'( A,A,T40,A,T60,1(1X,E20.12) )') ' INITIAL PRESSURE', & 00724 '[bar]', '= ', pv_scalar 00725 IF ( simpar%ensemble == nph_uniaxial_ensemble .OR.& 00726 simpar%ensemble == nph_uniaxial_damped_ensemble)& 00727 WRITE (iw,'( A,A,T40,A,T60,1(1X,E20.12) )') ' INITIAL HUGONIOT CONSTRAINT', & 00728 '[K]', '= ', hugoniot 00729 IF ( simpar%ensemble == nph_uniaxial_ensemble .OR.& 00730 simpar%ensemble == nph_uniaxial_damped_ensemble)& 00731 WRITE (iw,'( A,A,T40,A,T60,1(1X,E20.12) )') ' INITIAL E0', & 00732 '[hartree]', '= ', simpar%e0 00733 WRITE (iw,'( A,A,T40,A,T60,1(1X,E20.12) )') ' INITIAL VOLUME', & 00734 '[bohr^3]', '= ', cell%deth 00735 WRITE (iw,'( A,A,T29,A,T33,3(1X,E15.7) )') ' INITIAL CELL LNTHS', & 00736 '[bohr]', '= ', abc(1), abc(2), abc(3) 00737 WRITE (iw,'( A,A,T29,A,T33,3(1X,E15.7) )') ' INITIAL CELL ANGLS', & 00738 '[deg]', '= ', cell_angle(3), cell_angle(2), cell_angle(1) 00739 WRITE (iw,'( A,A )') ' ******************************** ', & 00740 'GO CP2K GO! **********************************' 00741 ELSE 00742 ! Write seuquential values of quantities of interest 00743 WRITE (iw,'(/,T2,A)') REPEAT('*',79) 00744 WRITE (iw,'(T2,A,T61,A20)')& 00745 'ENSEMBLE TYPE = ',ADJUSTR(TRIM(enum_i2c(enum,simpar%ensemble,error=error))) 00746 WRITE (iw,'(T2,A,T71,I10)')& 00747 'STEP NUMBER = ', itimes 00748 IF (simpar%variable_dt) THEN 00749 WRITE (iw,'(T2,A,T60,1X,F20.6)')& 00750 'TIME STEP [fs] = ',dt*femtoseconds 00751 END IF 00752 WRITE (iw,'(T2,A,T60,1X,F20.6)')& 00753 'TIME [fs] = ',time*femtoseconds 00754 WRITE (iw,'(T2,A,T60,1X,E20.12)')& 00755 'CONSERVED QUANTITY [hartree] = ',md_ener%constant 00756 WRITE (iw,'(/,T47,A,T73,A)') 'INSTANTANEOUS','AVERAGES' 00757 WRITE (iw,'(T2,A,T39,2(1X,F20.2))')& 00758 'CPU TIME [s] = ',used_time,averages%avecpu 00759 WRITE (iw,'(T2,A,T39,2(1X,E20.12))')& 00760 'ENERGY DRIFT PER ATOM [K] = ',econs,averages%econs 00761 WRITE (iw,'(T2,A,T39,2(1X,E20.12))')& 00762 'POTENTIAL ENERGY[hartree] = ',md_ener%epot,averages%avepot 00763 IF (simpar%ensemble /= reftraj_ensemble) THEN 00764 IF (.NOT.qmmm) THEN 00765 ! No QM/MM info 00766 WRITE (iw,'(T2,A,T39,2(1X,E20.12))')& 00767 'KINETIC ENERGY [hartree] = ',md_ener%ekin,averages%avekin 00768 WRITE (iw,'(T2,A,T39,2(1X,F20.3))')& 00769 'TEMPERATURE [K] = ',md_ener%temp_part,averages%avetemp 00770 ELSE 00771 WRITE (iw,'( A,A,T31,A,T39,2(1X,E20.12) )') ' TOTAL KINETIC ENERGY', & 00772 '[hartree]', '= ', md_ener%ekin, averages%avekin 00773 WRITE (iw,'( A,A,T31,A,T39,2(1X,E20.12) )') ' QM KINETIC ENERGY', & 00774 '[hartree]', '= ', md_ener%ekin_qm, averages%avekin_qm 00775 WRITE (iw,'( A,A,T31,A,T39,2(1X,F20.3) )') ' TOTAL TEMPERATURE', & 00776 '[K]', '= ', md_ener%temp_part, averages%avetemp 00777 WRITE (iw,'( A,A,T31,A,T39,2(1X,F20.3) )') ' QM TEMPERATURE', & 00778 '[K]', '= ', md_ener%temp_qm, averages%avetemp_qm 00779 END IF 00780 END IF 00781 IF (virial%pv_availability) THEN 00782 WRITE (iw,'(T2,A,T39,2(1X,E20.12))')& 00783 'PRESSURE [bar] = ',pv_scalar,averages%avepress 00784 END IF 00785 IF ( simpar%ensemble == nph_uniaxial_ensemble .OR.& 00786 simpar%ensemble == nph_uniaxial_damped_ensemble) THEN 00787 WRITE (iw,'( A,A,T31,A,T39,2(1X,E20.12) )') ' P_XX', & 00788 '[bar]', '= ', pv_xx, averages%avepxx 00789 WRITE (iw,'( A,A,T31,A,T39,2(1X,E20.12) )') ' HUGONIOT', & 00790 '[K]', '= ', hugoniot/3._dp/nat*kelvin,& 00791 averages%avehugoniot/3._dp/nat*kelvin 00792 END IF 00793 IF ( simpar%ensemble == nph_uniaxial_ensemble .OR.& 00794 simpar%ensemble == nph_uniaxial_damped_ensemble.OR.& 00795 simpar%ensemble == npt_i_ensemble.OR.& 00796 simpar%ensemble == npt_f_ensemble.OR.& 00797 simpar%ensemble == npe_i_ensemble.OR.& 00798 simpar%ensemble == npe_f_ensemble) THEN 00799 WRITE (iw,'( A,A,T31,A,T39,2(1X,E20.12) )') ' BAROSTAT TEMP', & 00800 '[K]', '= ', md_ener%temp_baro, averages%avetemp_baro 00801 WRITE (iw,'( A,A,T31,A,T39,2(1X,E20.12) )') ' VOLUME', & 00802 '[bohr^3]', '= ', cell%deth, averages%avevol 00803 WRITE (iw,'( A,A,T31,A,T33,3(1X,E15.7) )') ' CELL LNTHS', & 00804 '[bohr]', '= ', abc(1), abc(2), abc(3) 00805 WRITE (iw,'( A,A,T31,A,T33,3(1X,E15.7) )') ' AVE. CELL LNTHS', & 00806 '[bohr]', '= ', averages%aveca, averages%avecb, & 00807 averages%avecc 00808 END IF 00809 IF (simpar%ensemble==npt_f_ensemble .OR. & 00810 simpar%ensemble == npe_f_ensemble) THEN 00811 WRITE (iw,'( A,A,T31,A,T33,3(1X,E15.7) )') ' CELL ANGLS', & 00812 '[deg]', '= ', cell_angle(3), cell_angle(2), cell_angle(1) 00813 WRITE (iw,'( A,A,T31,A,T33,3(1X,E15.7) )') ' AVE. CELL ANGLS', & 00814 '[deg]', '= ', averages%aveal, averages%avebe, averages%avega 00815 END IF 00816 IF(simpar%ensemble==reftraj_ensemble) THEN 00817 IF(reftraj%info%msd) THEN 00818 IF (reftraj%msd%msd_kind) & 00819 WRITE(iw,'(/,A,T50,3f10.5)') ' COM displacement (dx,dy,dz) [angstrom]: ', & 00820 reftraj%msd%drcom(1:3)*angstrom 00821 END IF 00822 END IF 00823 WRITE (iw,'(T2,A,/)') REPEAT('*',79) 00824 END IF 00825 END IF 00826 CALL section_release(section, error) 00827 CALL cp_print_key_finished_output(iw,logger,motion_section,& 00828 "MD%PRINT%PROGRAM_RUN_INFO", error=error) 00829 END SUBROUTINE md_write_info_low 00830 00831 END MODULE md_energies
1.7.3