mirror of
https://github.com/cp2k/cp2k.git
synced 2026-07-28 14:15:19 -04:00
MD: implement temperature function (#5002)
This commit is contained in:
HE Zilong
GitHub
parent
a78533eb1f
commit
9f397624dc
4 changed files with 235 additions and 67 deletions
|
|
@ -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 <integer>", &
|
||||
CALL keyword_create(keyword, __LOCATION__, name="STEP_START_VAL", &
|
||||
description="The starting step value for the MD", usage="STEP_START_VAL <integer>", &
|
||||
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 <real>", default_r_val=cp_unit_to_cp2k(value=0.0_dp, unit_str="fs"), &
|
||||
usage="TIME_START_VAL <real>", 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 <real>", default_r_val=0.0_dp, &
|
||||
usage="ECONS_START_VAL <real>", 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)
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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'
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue