diff --git a/src/motion/input_cp2k_md.F b/src/motion/input_cp2k_md.F index 60abf1ac9a..959680db9a 100644 --- a/src/motion/input_cp2k_md.F +++ b/src/motion/input_cp2k_md.F @@ -22,6 +22,7 @@ MODULE input_cp2k_md low_print_level,& medium_print_level USE cp_units, ONLY: cp_unit_to_cp2k + USE fparser, ONLY: docf USE input_constants, ONLY: & isokin_ensemble, langevin_ensemble, md_init_default, md_init_vib, npe_f_ensemble, & npe_i_ensemble, nph_ensemble, nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, & @@ -40,14 +41,16 @@ MODULE input_cp2k_md section_create,& section_release,& section_type - USE input_val_types, ONLY: integer_t,& + USE input_val_types, ONLY: char_t,& + integer_t,& lchar_t,& real_t USE kinds, ONLY: dp USE reftraj_types, ONLY: REFTRAJ_EVAL_ENERGY,& REFTRAJ_EVAL_ENERGY_FORCES,& REFTRAJ_EVAL_NONE - USE string_utilities, ONLY: s2a + USE string_utilities, ONLY: newline,& + s2a #include "../base/base_uses.f90" IMPLICIT NONE @@ -106,99 +109,99 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="steps", & + CALL keyword_create(keyword, __LOCATION__, name="STEPS", & description="The number of MD steps to perform, counting from step_start_val. ", & - usage="steps 100", default_i_val=3) + usage="STEPS 100", default_i_val=3) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="max_steps", & + CALL keyword_create(keyword, __LOCATION__, name="MAX_STEPS", & description="The number of MD steps to perform, counting from step 1", & - usage="max_steps 100", default_i_val=1000000000) + usage="MAX_STEPS 100", default_i_val=1000000000) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="timestep", & + CALL keyword_create(keyword, __LOCATION__, name="TIMESTEP", & description="The length of an integration step (in case RESPA the large TIMESTEP)", & - usage="timestep 1.0", default_r_val=cp_unit_to_cp2k(value=0.5_dp, unit_str="fs"), & + usage="TIMESTEP 1.0", default_r_val=cp_unit_to_cp2k(value=0.5_dp, unit_str="fs"), & unit_str="fs") CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="step_start_val", & - description="The starting step value for the MD", usage="step_start_val ", & + CALL keyword_create(keyword, __LOCATION__, name="STEP_START_VAL", & + description="The starting step value for the MD", usage="STEP_START_VAL ", & default_i_val=0) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="time_start_val", & + CALL keyword_create(keyword, __LOCATION__, name="TIME_START_VAL", & description="The starting timer value for the MD", & - usage="time_start_val ", default_r_val=cp_unit_to_cp2k(value=0.0_dp, unit_str="fs"), & + usage="TIME_START_VAL ", default_r_val=cp_unit_to_cp2k(value=0.0_dp, unit_str="fs"), & unit_str="fs") CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="econs_start_val", & + CALL keyword_create(keyword, __LOCATION__, name="ECONS_START_VAL", & description="The starting value of the conserved quantity", & - usage="econs_start_val ", default_r_val=0.0_dp, & + usage="ECONS_START_VAL ", default_r_val=0.0_dp, & unit_str="hartree") CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="temperature", & + CALL keyword_create(keyword, __LOCATION__, name="TEMPERATURE", & description="The temperature in K used to initialize "// & "the velocities with init and pos restart, and in the NPT/NVT simulations", & - usage="temperature 325.0", default_r_val=cp_unit_to_cp2k(value=300.0_dp, unit_str="K"), & + usage="TEMPERATURE 325.0", default_r_val=cp_unit_to_cp2k(value=300.0_dp, unit_str="K"), & unit_str="K") CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="temp_tol", & + CALL keyword_create(keyword, __LOCATION__, name="TEMP_TOL", & variants=s2a("temp_to", "temperature_tolerance"), & description="The maximum accepted deviation of the (global) temperature "// & - "from the desired target temperature before a rescaling of the velocites "// & + "from the desired target temperature before a rescaling of the velocities "// & "is performed. If it is 0 no rescaling is performed. NOTE: This keyword is "// & "obsolescent; Using a CSVR thermostat with a short timeconstant is "// & "recommended as a better alternative.", & - usage="temp_tol 0.0", default_r_val=0.0_dp, unit_str='K') + usage="TEMP_TOL 0.0", default_r_val=0.0_dp, unit_str='K') CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="temp_kind", & + CALL keyword_create(keyword, __LOCATION__, name="TEMP_KIND", & description="Compute the temperature per each kind separately", & - usage="temp_kind LOGICAL", & + usage="TEMP_KIND LOGICAL", & default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="scale_temp_kind", & + CALL keyword_create(keyword, __LOCATION__, name="SCALE_TEMP_KIND", & description="When necessary rescale the temperature per each kind separately", & - usage="scale_temp_kind LOGICAL", & + usage="SCALE_TEMP_KIND LOGICAL", & default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="comvel_tol", & + CALL keyword_create(keyword, __LOCATION__, name="COMVEL_TOL", & description="The maximum accepted velocity of the center of mass. "// & "With Shell-Model, comvel may drift if MD%THERMOSTAT%REGION /= GLOBAL ", & - usage="comvel_tol 0.1", type_of_var=real_t, n_var=1, unit_str="bohr*au_t^-1") + usage="COMVEL_TOL 0.1", type_of_var=real_t, n_var=1, unit_str="bohr*au_t^-1") CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="angvel_tol", & + CALL keyword_create(keyword, __LOCATION__, name="ANGVEL_TOL", & description="The maximum accepted angular velocity. This option is ignored "// & "when the system is periodic. Removes the components of the velocities that "// & "project on the external rotational degrees of freedom.", & - usage="angvel_tol 0.1", type_of_var=real_t, n_var=1, unit_str="bohr*au_t^-1") + usage="ANGVEL_TOL 0.1", type_of_var=real_t, n_var=1, unit_str="bohr*au_t^-1") CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="angvel_zero", & + CALL keyword_create(keyword, __LOCATION__, name="ANGVEL_ZERO", & description="Set the initial angular velocity to zero. This option is ignored "// & "when the system is periodic or when initial velocities are defined. Technically, "// & "the part of the random initial velocities that projects on the external "// & "rotational degrees of freedom is subtracted.", & - usage="angvel_zero LOGICAL", & + usage="ANGVEL_ZERO LOGICAL", & default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) @@ -940,11 +943,12 @@ CONTAINS TYPE(section_type), POINTER :: section TYPE(keyword_type), POINTER :: keyword - TYPE(section_type), POINTER :: print_key, region_section, subsection + TYPE(section_type), POINTER :: print_key, region_section, subsection, & + tfunc_section CPASSERT(.NOT. ASSOCIATED(section)) - CALL section_create(section, __LOCATION__, name="thermal_region", & + CALL section_create(section, __LOCATION__, name="THERMAL_REGION", & description="Define regions where different initialization and control "// & "of the temperature is used. When MOTION%MD%ENSEMBLE is set to LANGEVIN, "// & "this section controls if the atoms defined inside and outside the "// & @@ -955,9 +959,9 @@ CONTAINS n_keywords=0, n_subsections=1, repeats=.FALSE.) NULLIFY (region_section) - NULLIFY (keyword, subsection) + NULLIFY (keyword, subsection, tfunc_section) - CALL keyword_create(keyword, __LOCATION__, name="force_rescaling", & + CALL keyword_create(keyword, __LOCATION__, name="FORCE_RESCALING", & description="Control the rescaling ot the velocities in all the regions, "// & "according to the temperature assigned to each reagion, when "// & "RESTART_VELOCITY in EXT_RESTART is active.", & @@ -965,18 +969,21 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="do_langevin_default", & + CALL keyword_create(keyword, __LOCATION__, name="DO_LANGEVIN_DEFAULT", & description="If ENSEMBLE is set to LANGEVIN, controls whether the "// & "atoms NOT defined in the thermal regions to undergo langevin MD "// & "or not. If not, then the atoms will undergo NVE Born-Oppenheimer MD.", & - usage="do_langevin_default .FALSE.", & + usage="DO_LANGEVIN_DEFAULT .FALSE.", & default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) CALL section_create(region_section, __LOCATION__, name="DEFINE_REGION", & - description="This section provides the possibility to define arbitrary region ", & - n_keywords=3, n_subsections=0, repeats=.TRUE.) + description="Define arbitrary thermal region with distinct "// & + "temperature controls. For thermostat regions, see section "// & + "MOTION%MD%THERMOSTAT%DEFINE_REGION. Typically the settings "// & + "should be the same.", & + n_keywords=3, n_subsections=1, repeats=.TRUE.) NULLIFY (keyword) CALL keyword_create(keyword, __LOCATION__, name="LIST", & @@ -986,43 +993,101 @@ CONTAINS CALL section_add_keyword(region_section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="temperature", & + CALL keyword_create(keyword, __LOCATION__, name="TEMPERATURE", & description="The temperature in K used to initialize the velocities "// & - "of the atoms in this region ", & - usage="temperature 5.0", & + "of the atoms in this region. It overrides the value evaluated by the "// & + "function defined in TFUNC section.", & + usage="TEMPERATURE 5.0", & default_r_val=cp_unit_to_cp2k(value=0.0_dp, unit_str="K"), & unit_str="K") CALL section_add_keyword(region_section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="temp_tol", & - description="Maximum accepted temperature deviation from the expected "// & - "value for this region. If temp_tol=0 no rescaling is performed", & - usage="temp_tol 0.0", & + CALL keyword_create(keyword, __LOCATION__, name="TEMP_TOL", & + description="Maximum accepted temperature deviation from the target "// & + "value for this region, which is defined by TEMPERATURE or evaluated "// & + "with the function defined in TFUNC section. Default value is 0 that "// & + "turns off rescaling.", & + usage="TEMP_TOL 0.0", & default_r_val=0.0_dp, unit_str='K') CALL section_add_keyword(region_section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="do_langevin", & + CALL keyword_create(keyword, __LOCATION__, name="DO_LANGEVIN", & description="When ENSEMBLE is set to LANGEVIN, Controls whether "// & "the atoms in the thermal region should undergo Langevin MD. If "// & "not, then they will undergo NVE Born-Oppenheimer MD.", & - usage="do_langevin .TRUE.", & + usage="DO_LANGEVIN .TRUE.", & default_l_val=.TRUE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(region_section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="noisy_gamma_region", & + CALL keyword_create(keyword, __LOCATION__, name="NOISY_GAMMA_REGION", & description="Special imaginary Langevin Friction term"// & " for Langevin Dynamics with noisy forces for the atoms in this region."// & " When set, overrides the general value set by NOISY_GAMMA in the MOTION%MD%LANGEVIN section."// & " When unset for a defined region, the general NOISY_GAMMA value applies.", & - citations=[Kuhne2007], usage="noisy_gamma_region 4.0E-5", & + citations=[Kuhne2007], usage="NOISY_GAMMA_REGION 4.0E-5", & type_of_var=real_t, & unit_str="fs^-1") CALL section_add_keyword(region_section, keyword) CALL keyword_release(keyword) + ! Define temperature function start + CALL section_create(tfunc_section, __LOCATION__, name="TFUNC", & + description="This section specifies the target temperature by using "// & + "the time step and the temperature of the last step, allowing for "// & + "simple region-wise simulated annealing (heating/cooling) controls. "// & + "When the absolute difference between target and actual temperature "// & + "is larger than TEMP_TOL value, the atomic velocities in this region "// & + "are rescaled to match the target temperature. However, note that "// & + "the actual temperature is subject to thermostat coupling as well "// & + "and a small restrictive TEMP_TOL together with a rapidly varying "// & + "target temperature function leaves little time for the system to "// & + "adjust itself, possibly showing unphysical relaxation behaviors. "// & + newline//"Further notes about temperature. "//docf(), & + n_keywords=4, n_subsections=0, repeats=.FALSE.) + + NULLIFY (keyword) + CALL keyword_create(keyword, __LOCATION__, name="FUNCTION", & + description="Specifies the functional form of the temperature "// & + "in Kelvin in mathematical notation. Two variables are available "// & + "with one denoted S for the integer time step and another denoted "// & + "T for the temperature of the last step in Kelvin respectively. When "// & + "DEFINE_REGION%TEMPERATURE is not explicitly set, the function "// & + "will be evaluated at S set in MD%STEP_START_VAL and T = 0.0 "// & + "to initialize the velocities of atoms in the region.", & + usage="FUNCTION A*T+(1-A)*B*(2/(1+EXP(-0.001*S))-1)", & + type_of_var=lchar_t, n_var=1) + CALL section_add_keyword(tfunc_section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="PARAMETERS", & + description="Defines the parameters of the functional form", & + usage="PARAMETERS A B", type_of_var=char_t, & + n_var=-1, repeats=.TRUE.) + CALL section_add_keyword(tfunc_section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="UNITS", & + description="Optionally, allows to define valid CP2K unit strings for each parameter value. "// & + "It is assumed that the corresponding parameter value is specified in this unit.", & + usage="UNITS ", type_of_var=char_t, & + n_var=-1, repeats=.TRUE.) + CALL section_add_keyword(tfunc_section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="VALUES", & + description="Defines the values of the parameter of the functional form", & + usage="VALUES 0.75 3.00E+2", type_of_var=real_t, & + n_var=-1, repeats=.TRUE., unit_str="internal_cp2k") + CALL section_add_keyword(tfunc_section, keyword) + CALL keyword_release(keyword) + + CALL section_add_subsection(region_section, tfunc_section) + CALL section_release(tfunc_section) + ! Define temperature function end + CALL section_add_subsection(section, region_section) CALL section_release(region_section) @@ -1032,7 +1097,13 @@ CONTAINS n_keywords=1, n_subsections=0, repeats=.FALSE.) CALL cp_print_key_section_create(print_key, __LOCATION__, "TEMPERATURE", & - description="Controls output of temperature per region.", & + description="Controls the output of the temperature per "// & + "region. The first three columns are: step number, time "// & + "(in fs), and the temperature of atoms not in any temperature "// & + "region (or zero if all atoms are in some temperature region). "// & + "Subsequent columns are, in pairs, the target value "// & + "(calculated by TFUNC or specified by TEMPERATURE) "// & + "and the actual value of the temperature for each region.", & print_level=high_print_level, common_iter_levels=1, & filename="") CALL section_add_subsection(subsection, print_key) diff --git a/src/motion/md_run.F b/src/motion/md_run.F index b57aaecb99..ff966f7ff3 100644 --- a/src/motion/md_run.F +++ b/src/motion/md_run.F @@ -87,7 +87,8 @@ MODULE md_run simpar_type USE thermal_region_types, ONLY: thermal_regions_type USE thermal_region_utils, ONLY: create_thermal_regions,& - print_thermal_regions_langevin + print_thermal_regions_langevin,& + update_thermal_regions_temp_expected USE thermostat_methods, ONLY: create_thermostats USE thermostat_types, ONLY: thermostats_type USE velocity_verlet_control, ONLY: velocity_verlet @@ -373,8 +374,6 @@ CONTAINS IF (ASSOCIATED(force_env%qs_env)) THEN CALL qs_env_time_update(force_env%qs_env, time, itimes) -!deb force_env%qs_env%sim_time = time -!deb force_env%qs_env%sim_step = itimes END IF ! Warm-up engines for metadynamics IF (ASSOCIATED(force_env%meta_env)) THEN @@ -456,8 +455,6 @@ CONTAINS !needed when electric field fields are applied IF (ASSOCIATED(force_env%qs_env)) THEN CALL qs_env_time_update(force_env%qs_env, time, itimes) -!deb force_env%qs_env%sim_time = time -!deb force_env%qs_env%sim_step = itimes END IF IF (ehrenfest_md) force_env%qs_env%rtp%istep = istep @@ -523,6 +520,10 @@ CONTAINS IF (simpar%ensemble /= reftraj_ensemble) THEN CALL md_energy(md_env, md_ener) + IF (simpar%do_thermal_region) THEN ! Optionally update expected temperature for each region + CALL get_md_env(md_env, thermal_regions=thermal_regions) + CALL update_thermal_regions_temp_expected(itimes, thermal_regions) + END IF CALL temperature_control(simpar, md_env, md_ener, force_env, logger) CALL comvel_control(md_ener, force_env, md_section, logger) CALL angvel_control(md_ener, force_env, md_section, logger) diff --git a/src/motion/thermal_region_types.F b/src/motion/thermal_region_types.F index e2c30a0926..d870dd8671 100644 --- a/src/motion/thermal_region_types.F +++ b/src/motion/thermal_region_types.F @@ -15,7 +15,9 @@ MODULE thermal_region_types USE input_section_types, ONLY: section_vals_type - USE kinds, ONLY: dp + USE kinds, ONLY: default_path_length,& + default_string_length,& + dp #include "../base/base_uses.f90" IMPLICIT NONE @@ -36,9 +38,12 @@ MODULE thermal_region_types END TYPE thermal_regions_type TYPE thermal_region_type + CHARACTER(len=default_path_length) :: temperature_function = "" + CHARACTER(len=default_string_length), POINTER, DIMENSION(:) :: temperature_function_parameters => NULL() INTEGER :: region_index = 0, npart = 0 INTEGER, DIMENSION(:), POINTER :: part_index => NULL() REAL(KIND=dp) :: ekin = 0.0_dp, noisy_gamma_region = 0.0_dp, temperature = 0.0_dp, temp_expected = 0.0_dp, temp_tol = 0.0_dp + REAL(KIND=dp), DIMENSION(:), POINTER :: temperature_function_values => NULL() END TYPE thermal_region_type CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'thermal_region_types' diff --git a/src/motion/thermal_region_utils.F b/src/motion/thermal_region_utils.F index d2f17c8413..c57335f77e 100644 --- a/src/motion/thermal_region_utils.F +++ b/src/motion/thermal_region_utils.F @@ -28,8 +28,16 @@ MODULE thermal_region_utils cp_print_key_unit_nr USE cp_subsys_types, ONLY: cp_subsys_get,& cp_subsys_type + USE cp_units, ONLY: cp_unit_from_cp2k,& + cp_unit_to_cp2k USE force_env_types, ONLY: force_env_get,& force_env_type + USE force_fields_util, ONLY: get_generic_info + USE fparser, ONLY: EvalErrType,& + evalf,& + finalizef,& + initf,& + parsef USE input_constants, ONLY: langevin_ensemble,& npt_f_ensemble,& npt_i_ensemble,& @@ -58,7 +66,8 @@ MODULE thermal_region_utils PRIVATE PUBLIC :: create_thermal_regions, & print_thermal_regions_temperature, & - print_thermal_regions_langevin + print_thermal_regions_langevin, & + update_thermal_regions_temp_expected CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'thermal_region_utils' @@ -81,18 +90,20 @@ CONTAINS TYPE(force_env_type), POINTER :: force_env CHARACTER(LEN=default_string_length) :: my_region - INTEGER :: i, il, ipart, ireg, nlist, nregions + INTEGER :: i, il, ipart, ireg, itimes, nlist, & + nregions INTEGER, DIMENSION(:), POINTER :: tmplist LOGICAL :: apply_thermostat, do_langevin, & do_langevin_default, do_read_ngr, & - explicit + explicit, use_t, use_tfunc REAL(KIND=dp) :: temp, temp_tol TYPE(cp_subsys_type), POINTER :: subsys TYPE(particle_list_type), POINTER :: particles - TYPE(section_vals_type), POINTER :: region_sections, thermal_region_section + TYPE(section_vals_type), POINTER :: region_sections, tfunc_section, & + thermal_region_section TYPE(thermal_region_type), POINTER :: t_region - NULLIFY (region_sections, t_region, thermal_region_section, particles, subsys, tmplist) + NULLIFY (region_sections, t_region, tfunc_section, thermal_region_section, particles, subsys, tmplist) ALLOCATE (thermal_regions) CALL allocate_thermal_regions(thermal_regions) thermal_region_section => section_vals_get_subs_vals(md_section, "THERMAL_REGION") @@ -135,8 +146,11 @@ CONTAINS END IF DO ireg = 1, nregions NULLIFY (t_region) + CALL integer_to_string(ireg, my_region) t_region => thermal_regions%thermal_region(ireg) t_region%region_index = ireg + + ! Determine t_region%npart CALL section_vals_val_get(region_sections, "LIST", & i_rep_section=ireg, n_rep_val=nlist) NULLIFY (t_region%part_index) @@ -160,9 +174,39 @@ CONTAINS END IF END DO END DO - CALL section_vals_val_get(region_sections, "TEMPERATURE", i_rep_section=ireg, & - r_val=temp) - t_region%temp_expected = temp + + ! Determine t_region%temp_expected + tfunc_section => section_vals_get_subs_vals(region_sections, "TFUNC", i_rep_section=ireg) + CALL section_vals_val_get(region_sections, "TEMPERATURE", i_rep_section=ireg, explicit=use_t) + IF (use_t) CALL section_vals_val_get(region_sections, "TEMPERATURE", i_rep_section=ireg, r_val=temp) + CALL section_vals_get(tfunc_section, explicit=use_tfunc) + IF (use_tfunc) THEN + CALL get_generic_info(tfunc_section, "FUNCTION", & + t_region%temperature_function, & + t_region%temperature_function_parameters, & + t_region%temperature_function_values, & + input_variables=["S", "T"], i_rep_sec=1) + CALL initf(1) + CALL parsef(1, TRIM(t_region%temperature_function), & + t_region%temperature_function_parameters) + IF (.NOT. use_t) THEN ! Evaluate function for initial temperature + CALL section_vals_val_get(md_section, "STEP_START_VAL", i_val=itimes) + t_region%temperature_function_values(1) = itimes*1.0_dp + t_region%temperature_function_values(2) = 0.0_dp + temp = cp_unit_to_cp2k(evalf(1, t_region%temperature_function_values), "K") + CPASSERT(EvalErrType <= 0) + END IF + CALL finalizef() + END IF + IF (temp >= 0.0_dp) THEN ! Check for negative temperature in Kelvin + t_region%temp_expected = temp + ELSE + CALL cp_abort(__LOCATION__, & + "The initial temperature for region "//TRIM(my_region)//" is "// & + TRIM(ADJUSTL(cp_to_string(cp_unit_from_cp2k(temp, "K"))))// & + " K, which is negative and unphysical!") + END IF + CALL section_vals_val_get(region_sections, "TEMP_TOL", i_rep_section=ireg, & r_val=temp_tol) t_region%temp_tol = temp_tol @@ -172,7 +216,6 @@ CONTAINS r_val=t_region%noisy_gamma_region) IF (simpar%ensemble == langevin_ensemble) THEN IF (.NOT. do_langevin) THEN - CALL integer_to_string(ireg, my_region) CALL cp_warn(__LOCATION__, & "You provided NOISY_GAMMA_REGION but atoms in thermal region "//TRIM(my_region)// & " will not undergo Langevin MD. "// & @@ -201,6 +244,50 @@ CONTAINS END SUBROUTINE create_thermal_regions +! ************************************************************************************************** +!> \brief update_thermal_regions_temp_expected +!> \param itimes : iteration number of the time step +!> \param thermal_regions : thermal regions type contains information +!> about the regions +!> \date 03.2026 +!> \par Note on units: t_region%temperature is updated by +!> SUBROUTINE scale_velocity_region in md_vel_utils.F +!> which is in Kelvin, but t_region%temp_expected is +!> in cp2k internal units. See "Note the unit conversion +!> here" in print_thermal_regions_temperature below. +!> \author HE Zilong +! ************************************************************************************************** + SUBROUTINE update_thermal_regions_temp_expected(itimes, thermal_regions) + INTEGER, INTENT(IN) :: itimes + TYPE(thermal_regions_type), POINTER :: thermal_regions + + CHARACTER(LEN=default_string_length) :: my_region + INTEGER :: ireg + TYPE(thermal_region_type), POINTER :: t_region + + DO ireg = 1, thermal_regions%nregions + NULLIFY (t_region) + t_region => thermal_regions%thermal_region(ireg) + IF (LEN_TRIM(t_region%temperature_function) > 0) THEN + CALL initf(1) + CALL parsef(1, TRIM(t_region%temperature_function), t_region%temperature_function_parameters) + t_region%temperature_function_values(1) = itimes*1.0_dp + t_region%temperature_function_values(2) = t_region%temperature + t_region%temp_expected = cp_unit_to_cp2k(evalf(1, t_region%temperature_function_values), "K") + CPASSERT(EvalErrType <= 0) + CALL finalizef() + IF (t_region%temp_expected < 0.0_dp) THEN + CALL integer_to_string(ireg, my_region) + CALL cp_abort(__LOCATION__, & + "The temperature evaluated for region "//TRIM(my_region)//" is "// & + TRIM(ADJUSTL(cp_to_string(cp_unit_from_cp2k(t_region%temp_expected, "K"))))// & + " K, which is negative and unphysical!") + END IF + END IF + END DO + + END SUBROUTINE update_thermal_regions_temp_expected + ! ************************************************************************************************** !> \brief print_thermal_regions_temperature !> \param thermal_regions : thermal regions type contains information @@ -224,7 +311,7 @@ CONTAINS CHARACTER(LEN=default_string_length) :: fmd INTEGER :: ireg, nregions, unit LOGICAL :: new_file - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: temp + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: temp, temp_tfunc TYPE(cp_logger_type), POINTER :: logger TYPE(section_vals_type), POINTER :: print_key @@ -244,15 +331,19 @@ CONTAINS END IF nregions = thermal_regions%nregions ALLOCATE (temp(0:nregions)) + ALLOCATE (temp_tfunc(1:nregions)) temp = 0.0_dp temp(0) = thermal_regions%temp_reg0 DO ireg = 1, nregions + ! Note the unit conversion here temp(ireg) = thermal_regions%thermal_region(ireg)%temperature + temp_tfunc(ireg) = cp_unit_from_cp2k(thermal_regions%thermal_region(ireg)%temp_expected, "K") END DO - fmd = "(I10,F20.3,"//TRIM(ADJUSTL(cp_to_string(nregions + 1)))//"F20.6)" + fmd = "(I10,F20.3,"//TRIM(ADJUSTL(cp_to_string(2*nregions + 1)))//"F20.8)" fmd = TRIM(fmd) - WRITE (UNIT=unit, FMT=fmd) itimes, time, temp(0:nregions) + WRITE (UNIT=unit, FMT=fmd) itimes, time, temp(0), (temp_tfunc(ireg), temp(ireg), ireg=1, nregions) DEALLOCATE (temp) + DEALLOCATE (temp_tfunc) END IF CALL cp_print_key_finished_output(unit, logger, thermal_regions%section, "PRINT%TEMPERATURE") END IF