|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00012 MODULE tamc_run 00013 !!!!!!!!!! 00014 !!!!!!!!!! integrators 00015 USE atomic_kind_list_types, ONLY: atomic_kind_list_type 00016 USE atomic_kind_types, ONLY: atomic_kind_type 00017 USE averages_types, ONLY: average_quantities_type 00018 USE barostat_types, ONLY: barostat_type,& 00019 create_barostat_type,& 00020 release_barostat_type 00021 USE bibliography, ONLY: VandenCic2006 00022 USE cell_types, ONLY: cell_type 00023 USE colvar_methods, ONLY: colvar_eval_glob_f 00024 USE colvar_types, ONLY: HBP_colvar_id,& 00025 WC_colvar_id,& 00026 colvar_p_type 00027 USE constraint_fxd, ONLY: fix_atom_control 00028 USE cp_external_control, ONLY: external_control 00029 USE cp_output_handling, ONLY: cp_add_iter_level,& 00030 cp_iterate,& 00031 cp_p_file,& 00032 cp_print_key_finished_output,& 00033 cp_print_key_should_output,& 00034 cp_print_key_unit_nr,& 00035 cp_rm_iter_level 00036 USE cp_para_types, ONLY: cp_para_env_type 00037 USE cp_subsys_types, ONLY: cp_subsys_get,& 00038 cp_subsys_type 00039 USE cp_units, ONLY: cp_unit_from_cp2k 00040 USE distribution_1d_types, ONLY: distribution_1d_type 00041 USE f77_blas 00042 USE force_env_methods, ONLY: force_env_calc_energy_force 00043 USE force_env_types, ONLY: force_env_get,& 00044 force_env_type 00045 USE free_energy_types, ONLY: fe_env_create,& 00046 free_energy_type 00047 USE global_types, ONLY: global_environment_type 00048 USE input_constants, ONLY: & 00049 ehrenfest, langevin_ensemble, npe_f_ensemble, npe_i_ensemble, & 00050 nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, & 00051 npt_i_ensemble, reftraj_ensemble 00052 USE input_cp2k_check, ONLY: remove_restart_info 00053 USE input_cp2k_restarts, ONLY: write_restart 00054 USE input_section_types, ONLY: section_vals_get,& 00055 section_vals_get_subs_vals,& 00056 section_vals_remove_values,& 00057 section_vals_type,& 00058 section_vals_val_get,& 00059 section_vals_val_set 00060 USE kinds, ONLY: dp 00061 USE machine, ONLY: m_walltime 00062 USE mc_environment_types, ONLY: get_mc_env,& 00063 mc_env_create,& 00064 mc_env_release,& 00065 mc_environment_type,& 00066 set_mc_env 00067 USE mc_misc, ONLY: mc_averages_create,& 00068 mc_averages_release 00069 USE mc_move_control, ONLY: init_mc_moves,& 00070 mc_moves_release 00071 USE mc_types, ONLY: get_mc_par,& 00072 mc_averages_type,& 00073 mc_ekin_type,& 00074 mc_moves_type,& 00075 mc_simpar_type,& 00076 set_mc_par 00077 USE md_ener_types, ONLY: create_md_ener,& 00078 md_ener_type,& 00079 release_md_ener 00080 USE md_energies, ONLY: initialize_md_ener,& 00081 md_energy 00082 USE md_environment_types, ONLY: get_md_env,& 00083 md_env_create,& 00084 md_env_release,& 00085 md_environment_type,& 00086 set_md_env 00087 USE md_run, ONLY: qs_mol_dyn 00088 USE metadynamics_types, ONLY: meta_env_type,& 00089 metavar_type,& 00090 set_meta_env 00091 USE mol_kind_new_list_types, ONLY: mol_kind_new_list_type 00092 USE mol_new_list_types, ONLY: mol_new_list_type 00093 USE molecule_kind_types, ONLY: molecule_kind_type 00094 USE molecule_types_new, ONLY: global_constraint_type,& 00095 molecule_type 00096 USE parallel_rng_types, ONLY: UNIFORM,& 00097 create_rng_stream,& 00098 delete_rng_stream,& 00099 next_random_number,& 00100 rng_stream_type 00101 USE particle_list_types, ONLY: particle_list_type 00102 USE particle_types, ONLY: particle_type 00103 USE physcon, ONLY: boltzmann,& 00104 femtoseconds,& 00105 joule,& 00106 kelvin 00107 USE qmmm_util, ONLY: apply_qmmm_walls_reflective 00108 USE qs_environment_types, ONLY: get_qs_env 00109 USE qs_scf_post_gpw, ONLY: scf_post_calculation_gpw 00110 USE qs_scf_types, ONLY: qs_scf_env_type 00111 USE reference_manager, ONLY: cite_reference 00112 USE reftraj_types, ONLY: create_reftraj,& 00113 reftraj_type,& 00114 release_reftraj 00115 USE reftraj_util, ONLY: initialize_reftraj 00116 USE rt_propagation, ONLY: rt_prop_setup,& 00117 rt_write_input_restart 00118 USE rt_propagation_output, ONLY: rt_prop_output 00119 USE simpar_methods, ONLY: read_md_section 00120 USE simpar_types, ONLY: create_simpar_type,& 00121 release_simpar_type,& 00122 simpar_type 00123 USE string_utilities, ONLY: str_comp 00124 USE thermal_region_types, ONLY: release_thermal_regions,& 00125 thermal_regions_type 00126 USE thermal_region_utils, ONLY: create_thermal_regions 00127 USE thermostat_methods, ONLY: create_thermostats 00128 USE thermostat_types, ONLY: release_thermostats,& 00129 thermostats_type 00130 USE timings, ONLY: timeset,& 00131 timestop 00132 USE virial_methods, ONLY: virial_evaluate 00133 USE virial_types, ONLY: virial_type 00134 USE wiener_process, ONLY: create_wiener_process,& 00135 create_wiener_process_cv 00136 !!!!! monte carlo part 00137 #include "cp_common_uses.h" 00138 00139 IMPLICIT NONE 00140 00141 00142 PRIVATE 00143 00144 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tamc_run' 00145 00146 PUBLIC :: qs_tamc 00147 00148 CONTAINS 00149 00150 ! ***************************************************************************** 00155 SUBROUTINE qs_tamc(force_env, globenv,input_file_name,averages,error) 00156 00157 TYPE(force_env_type), POINTER :: force_env 00158 TYPE(global_environment_type), POINTER :: globenv 00159 CHARACTER(LEN=*), OPTIONAL :: input_file_name 00160 TYPE(average_quantities_type), 00161 OPTIONAL, POINTER :: averages 00162 TYPE(cp_error_type), INTENT(inout) :: error 00163 00164 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_tamc', 00165 routineP = moduleN//':'//routineN 00166 00167 CHARACTER(LEN=20) :: ensemble 00168 INTEGER :: handle, i, initialStep, iprint, isos, istep, j, md_stride, 00169 nmccycles, output_unit, rand2skip, run_type_id 00170 INTEGER, POINTER :: itimes 00171 LOGICAL :: check, ehrenfest_md, 00172 explicit, failure, ionode, 00173 my_rm_restart_info, save_mem, 00174 should_stop 00175 REAL(KIND=dp) :: auxRandom, inittime, rval, 00176 temp, time_iter_start, 00177 time_iter_stop 00178 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: An, fz, xieta, zbuff 00179 REAL(KIND=dp), ALLOCATABLE, 00180 DIMENSION(:, :) :: r 00181 REAL(KIND=dp), POINTER :: constant, time, used_time 00182 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 00183 TYPE(barostat_type), POINTER :: barostat 00184 TYPE(cell_type), POINTER :: cell, cell_ref 00185 TYPE(cp_logger_type), POINTER :: logger 00186 TYPE(cp_para_env_type), POINTER :: para_env 00187 TYPE(cp_subsys_type), POINTER :: subsys 00188 TYPE(distribution_1d_type), POINTER :: local_particles 00189 TYPE(free_energy_type), POINTER :: fe_env 00190 TYPE(mc_averages_type), POINTER :: MCaverages 00191 TYPE(mc_environment_type), POINTER :: mc_env 00192 TYPE(mc_moves_type), POINTER :: gmoves, moves 00193 TYPE(mc_simpar_type), POINTER :: mc_par 00194 TYPE(md_ener_type), POINTER :: md_ener 00195 TYPE(md_environment_type), POINTER :: md_env 00196 TYPE(meta_env_type), POINTER :: meta_env_saved 00197 TYPE(particle_list_type), POINTER :: particles 00198 TYPE(reftraj_type), POINTER :: reftraj 00199 TYPE(rng_stream_type), POINTER :: rng_stream, rng_stream_mc 00200 TYPE(section_vals_type), POINTER :: constraint_section, 00201 force_env_section, free_energy_section, fs_section, global_section, 00202 mc_section, md_section, motion_section, reftraj_section, 00203 subsys_section, work_section 00204 TYPE(simpar_type), POINTER :: simpar 00205 TYPE(thermal_regions_type), POINTER :: thermal_regions 00206 TYPE(thermostats_type), POINTER :: thermostats 00207 00208 initialStep=0 00209 inittime=0.0_dp 00210 00211 CALL timeset(routineN,handle) 00212 failure = .FALSE. 00213 my_rm_restart_info = .TRUE. 00214 NULLIFY(md_env, para_env,fs_section) 00215 para_env => force_env%para_env 00216 motion_section => section_vals_get_subs_vals(force_env%root_section,"MOTION",error=error) 00217 md_section => section_vals_get_subs_vals(motion_section,"MD",error=error) 00218 00219 ! Real call to MD driver - Low Level 00220 CALL md_env_create(md_env, md_section, para_env, force_env=force_env, error=error) 00221 IF (PRESENT(averages)) CALL set_md_env(md_env, averages=averages, error=error) 00222 00223 00224 00225 CPPrecondition(ASSOCIATED(globenv),cp_failure_level,routineP,error,failure) 00226 CPPrecondition(ASSOCIATED(force_env),cp_failure_level,routineP,error,failure) 00227 00228 failure=.FALSE. 00229 NULLIFY (particles, cell, cell_ref, simpar, itimes, used_time, subsys, & 00230 md_ener, thermostats, barostat, reftraj, force_env_section, & 00231 reftraj_section, work_section, atomic_kinds, & 00232 local_particles, time, fe_env, free_energy_section, & 00233 constraint_section, thermal_regions) 00234 logger => cp_error_get_logger(error) 00235 para_env => force_env%para_env 00236 00237 global_section => section_vals_get_subs_vals(force_env%root_section,"GLOBAL",error=error) 00238 free_energy_section =>section_vals_get_subs_vals(motion_section,"FREE_ENERGY",error=error) 00239 constraint_section =>section_vals_get_subs_vals(motion_section,"CONSTRAINT",error=error) 00240 CALL section_vals_val_get(global_section,"SAVE_MEM",l_val=save_mem,error=error) 00241 00242 CALL section_vals_val_get(global_section,"RUN_TYPE", i_val=run_type_id,error=error) 00243 IF(run_type_id==ehrenfest) CALL set_md_env(md_env, ehrenfest_md=.TRUE., error=error) 00244 00245 CALL create_simpar_type(simpar, error) 00246 force_env_section => force_env%force_env_section 00247 subsys_section => section_vals_get_subs_vals(force_env_section,"SUBSYS",error=error) 00248 CALL cp_add_iter_level(logger%iter_info,"MD",error=error) 00249 CALL cp_iterate(logger%iter_info,iter_nr=initialStep,error=error) 00250 ! Read MD section 00251 CALL read_md_section(simpar, motion_section, md_section, error) 00252 ! Setup print_keys 00253 simpar%info_constraint = cp_print_key_unit_nr(logger,constraint_section,& 00254 "CONSTRAINT_INFO",extension=".shakeLog",log_filename=.FALSE.,error=error) 00255 simpar%lagrange_multipliers = cp_print_key_unit_nr(logger,constraint_section,& 00256 "LAGRANGE_MULTIPLIERS",extension=".LagrangeMultLog",log_filename=.FALSE.,error=error) 00257 simpar%dump_lm = BTEST(cp_print_key_should_output(logger%iter_info,constraint_section,& 00258 "LAGRANGE_MULTIPLIERS",error=error),cp_p_file) 00259 00260 ! Create the structure for the md energies 00261 CALL create_md_ener(md_ener, error=error) 00262 CALL set_md_env(md_env, md_ener=md_ener, error=error) 00263 CALL release_md_ener(md_ener, error=error) 00264 00265 ! If requested setup Thermostats 00266 CALL create_thermostats(thermostats, md_section, force_env, simpar, para_env,& 00267 globenv, global_section, error ) 00268 00269 ! If requested setup Barostat 00270 CALL create_barostat_type(barostat, md_section, force_env, simpar, globenv, error ) 00271 00272 ! If requested setup different thermal regions 00273 CALL create_thermal_regions(thermal_regions, md_section, simpar, force_env, error ) 00274 00275 CALL set_md_env(md_env, thermostats=thermostats, barostat=barostat, thermal_regions=thermal_regions,error=error) 00276 00277 CALL get_md_env(md_env, ehrenfest_md=ehrenfest_md, error=error) 00278 00279 !If requested set up the REFTRAJ run 00280 IF(simpar%ensemble == reftraj_ensemble .AND. ehrenfest_md)& 00281 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 00282 "Ehrenfest MD does not support reftraj ensemble "//& 00283 CPSourceFileRef,& 00284 error,failure) 00285 IF(simpar%ensemble == reftraj_ensemble) THEN 00286 reftraj_section => section_vals_get_subs_vals(md_section,"REFTRAJ",error=error) 00287 CALL create_reftraj(reftraj, reftraj_section, para_env, error=error) 00288 CALL set_md_env(md_env, reftraj=reftraj, error=error) 00289 CALL release_reftraj(reftraj,error=error) 00290 END IF 00291 00292 CALL force_env_get(force_env, subsys=subsys, cell=cell, cell_ref=cell_ref,& 00293 force_env_section=force_env_section, error=error ) 00294 00295 00296 00297 00298 00299 ! Set V0 if needed 00300 IF (simpar%ensemble == nph_uniaxial_ensemble.OR.simpar%ensemble == nph_uniaxial_damped_ensemble) THEN 00301 CPPrecondition(ASSOCIATED(cell_ref),cp_failure_level,routineP,error,failure) 00302 IF ( simpar%v0 == 0._dp ) simpar%v0 = cell_ref%deth 00303 ENDIF 00304 00305 ! Initialize velocities possibly applying constraints at the zeroth MD step 00306 ! ! ! CALL section_vals_val_get(motion_section,"PRINT%RESTART%SPLIT_RESTART_FILE",& 00307 ! ! ! l_val=write_binary_restart_file,error=error) 00308 !! let us see if this created all the trouble 00309 ! CALL setup_velocities(force_env,simpar,globenv,md_env,md_section,constraint_section, & 00310 ! write_binary_restart_file,error) 00311 00312 ! Setup Free Energy Calculation (if required) 00313 CALL fe_env_create (fe_env, free_energy_section, error) 00314 CALL set_md_env(md_env=md_env, simpar=simpar, fe_env=fe_env, cell=cell,& 00315 force_env=force_env, error=error) 00316 00317 ! Possibly initialize Wiener processes 00318 IF (simpar%ensemble == langevin_ensemble) CALL create_wiener_process(md_env,error) 00319 time_iter_start=m_walltime() 00320 00321 CALL get_md_env(md_env, force_env=force_env, itimes=itimes, constant=constant,& 00322 md_ener=md_ener, t=time, used_time=used_time, error=error) 00323 00324 ! Attach the time counter of the meta_env to the one of the MD 00325 CALL set_meta_env(force_env%meta_env, time=time, error=error) 00326 ! Initialize the md_ener structure 00327 00328 force_env%meta_env%dt=force_env%meta_env%zdt 00329 CALL initialize_md_ener(md_ener, force_env, simpar, error=error) 00330 ! force_env%meta_env%dt=force_env%meta_env%zdt 00331 00332 00333 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! MC setup up 00334 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00335 00336 NULLIFY(mc_env,mc_par,rng_stream_mc,MCaverages) 00337 00338 CALL section_vals_get(force_env_section,n_repetition=isos,error=error) 00339 CPPostconditionNoFail(isos==1,cp_warning_level,routineP,error) 00340 ! set some values...will use get_globenv if that ever comes around 00341 00342 ! initialize the random numbers 00343 ! IF (para_env%ionode) THEN 00344 CALL create_rng_stream(rng_stream=rng_stream_mc,& 00345 name="Random numbers for monte carlo acc/rej",& 00346 distribution_type=UNIFORM,error=error) 00347 ! ENDIF 00348 !!!!! this shoudl go in a routine hmc_read 00349 00350 NULLIFY(mc_section) 00351 ALLOCATE(mc_par, stat=isos) ! do not forget to clean this mess 00352 CPPostcondition(isos==0,cp_failure_level,routineP,error,failure) 00353 00354 mc_section => section_vals_get_subs_vals(force_env%root_section,& 00355 "MOTION%MC",error=error) 00356 CALL section_vals_val_get(mc_section,"ENSEMBLE",& 00357 c_val=ensemble,error=error) 00358 CPPostcondition(str_comp(ensemble,"TRADITIONAL"),cp_failure_level,routineP,error,failure) 00359 CALL section_vals_val_get(mc_section,"NSTEP",& 00360 i_val=nmccycles,error=error) 00361 CPPostcondition(nmccycles>0,cp_failure_level,routineP,error,failure) 00362 CALL section_vals_val_get(mc_section,"IPRINT",& 00363 i_val=iprint,error=error) 00364 CALL section_vals_val_get(mc_section,"RANDOMTOSKIP",i_val=rand2skip,error=error) 00365 CPPostcondition(rand2skip>=0,cp_failure_level,routineP,error,failure) 00366 temp=cp_unit_from_cp2k(simpar%temp_ext,"K",error=error) 00367 ! 00368 00369 CALL set_mc_par(mc_par, ensemble=ensemble, nstep=nmccycles, iprint=iprint, temperature=temp, & 00370 beta=1.0_dp / temp / boltzmann * joule,exp_max_val = 0.9_dp*LOG(HUGE(0.0_dp)),& 00371 exp_min_val = 0.9_dp*LOG(TINY(0.0_dp)),max_val=HUGE(0.0_dp),min_val= 0.0_dp, & 00372 source=para_env%source, group=para_env%group, ionode=para_env%ionode, rand2skip=rand2skip) 00373 00374 output_unit = cp_logger_get_default_io_unit(logger) 00375 IF (output_unit > 0) THEN 00376 WRITE(output_unit,'(a,a)')"HMC| Hybrid Monte Carlo Scheme " 00377 WRITE(output_unit,'(a,a)')"HMC| Ensemble ", ADJUSTL(ensemble) 00378 WRITE(output_unit,'(a,i0)')"HMC| MC Cycles ", nmccycles 00379 WRITE(output_unit,'(a,i0,a)')"HMC| Print every ", iprint, " cycles" 00380 WRITE(output_unit,'(a,i0)')"HMC| Number of random numbers to skip ", rand2skip 00381 WRITE(output_unit,'(a,f16.8,a)')"HMC| Temperature ",temp, "K" 00382 ENDIF 00383 00384 CALL force_env_get( force_env, subsys=subsys, error=error ) 00385 00386 CALL cp_subsys_get(subsys,atomic_kinds=atomic_kinds,local_particles=local_particles,& 00387 particles=particles,error=error) 00388 00389 DO i=1,rand2skip 00390 auxRandom=next_random_number(rng_stream_mc,error=error) 00391 DO j=1,3*SIZE(particles%els) 00392 auxRandom=next_random_number(globenv%gaussian_rng_stream,error=error) 00393 ENDDO 00394 ENDDO 00395 00396 00397 00398 CALL mc_env_create ( mc_env, error = error ) 00399 CALL set_mc_env( mc_env,mc_par=mc_par,force_env=force_env) 00400 !!!!!!!end mc setup 00401 00402 ! Check for ensembles requiring the stress tensor - takes into account the possibility for 00403 ! multiple force_evals 00404 IF ( (simpar%ensemble==npt_i_ensemble).OR.& 00405 (simpar%ensemble==npt_f_ensemble).OR.& 00406 (simpar%ensemble==npe_f_ensemble).OR.& 00407 (simpar%ensemble==npe_i_ensemble).OR.& 00408 (simpar%ensemble==nph_uniaxial_ensemble).OR.& 00409 (simpar%ensemble==nph_uniaxial_damped_ensemble)) THEN 00410 check = force_env%virial%pv_availability 00411 CALL cp_assert(check,cp_failure_level,cp_assertion_failed,& 00412 routineP,"Virial evaluation not requested for this run in the input file! "//& 00413 "You may consider to switch on the virial evaluation with the keyword: STRESS_TENSOR."//& 00414 "Be sure the method you are using can compute the virial! "//& 00415 CPSourceFileRef,& 00416 error,failure) 00417 IF (ASSOCIATED(force_env%sub_force_env)) THEN 00418 DO i = 1, SIZE(force_env%sub_force_env) 00419 IF (ASSOCIATED(force_env%sub_force_env(i)%force_env)) THEN 00420 check = check .AND. force_env%sub_force_env(i)%force_env%virial%pv_availability 00421 END IF 00422 END DO 00423 END IF 00424 CALL cp_assert(check,cp_failure_level,cp_assertion_failed,& 00425 routineP,"Virial evaluation not requested for all the force_eval sections present in"//& 00426 " the input file! You have to switch on the virial evaluation with the keyword: STRESS_TENSOR "//& 00427 " in each force_eval section. Be sure the method you are using can compute the virial!"//& 00428 CPSourceFileRef,& 00429 error,failure) 00430 END IF 00431 00432 ! Computing Forces at zero MD step 00433 IF (simpar%ensemble /= reftraj_ensemble) THEN 00434 CALL section_vals_val_get(md_section,"STEP_START_VAL",i_val=itimes,error=error) 00435 CALL section_vals_val_get(md_section,"TIME_START_VAL",r_val=time,error=error) 00436 CALL section_vals_val_get(md_section,"ECONS_START_VAL",r_val=constant,error=error) 00437 CALL section_vals_val_set(md_section,"STEP_START_VAL",i_val=initialStep,error=error) 00438 CALL section_vals_val_set(md_section,"TIME_START_VAL",r_val=inittime,error=error) 00439 initialStep=itimes 00440 CALL cp_iterate(logger%iter_info,iter_nr=itimes,error=error) 00441 IF(save_mem) THEN 00442 work_section => section_vals_get_subs_vals(subsys_section,"VELOCITY",error=error) 00443 CALL section_vals_remove_values(work_section, error) 00444 work_section => section_vals_get_subs_vals(subsys_section,"SHELL_VELOCITY",error=error) 00445 CALL section_vals_remove_values(work_section, error) 00446 work_section => section_vals_get_subs_vals(subsys_section,"CORE_VELOCITY",error=error) 00447 CALL section_vals_remove_values(work_section, error) 00448 END IF 00449 00450 IF(ehrenfest_md)THEN 00451 CALL rt_prop_setup(force_env,error) 00452 force_env%qs_env%rtp%dt=simpar%dt 00453 ELSE 00454 ! CALL force_env_calc_energy_force (force_env, calc_force=.TRUE., error=error) 00455 meta_env_saved=> force_env%meta_env 00456 NULLIFY(force_env%meta_env) 00457 CALL force_env_calc_energy_force (force_env, calc_force=.FALSE., error=error) 00458 force_env%meta_env=>meta_env_saved 00459 END IF 00460 00461 IF(ASSOCIATED(force_env%qs_env))THEN 00462 ! force_env%qs_env%sim_time=time 00463 ! force_env%qs_env%sim_step=itimes 00464 force_env%qs_env%sim_time=0.0_dp 00465 force_env%qs_env%sim_step=0 00466 END IF 00467 ! Warm-up engines for metadynamics 00468 IF (ASSOCIATED(force_env%meta_env)) THEN 00469 IF(force_env%meta_env%langevin) THEN 00470 CALL create_wiener_process_cv(force_env%meta_env, error=error) 00471 NULLIFY(rng_stream) 00472 DO j=1, (rand2skip-1)/nmccycles 00473 DO i = 1 , force_env%meta_env%n_colvar 00474 rng_stream => force_env%meta_env%rng(i)%stream 00475 auxRandom=next_random_number(rng_stream,error=error) 00476 auxRandom=next_random_number(rng_stream,error=error) 00477 ENDDO 00478 ENDDO 00479 ENDIF 00480 ! IF (force_env%meta_env%well_tempered) THEN 00481 ! force_env%meta_env%wttemperature = simpar%temp_ext 00482 ! IF (force_env%meta_env%wtgamma>EPSILON(1._dp)) THEN 00483 ! dummy=force_env%meta_env%wttemperature*(force_env%meta_env%wtgamma-1._dp) 00484 ! IF (force_env%meta_env%delta_t>EPSILON(1._dp)) THEN 00485 ! check=ABS(force_env%meta_env%delta_t-dummy)<1.E+3_dp*EPSILON(1._dp) 00486 ! CALL cp_assert(check,cp_failure_level,cp_assertion_failed,routineP,& 00487 ! "Inconsistency between DELTA_T and WTGAMMA (both specified):"//& 00488 ! " please, verify that DELTA_T=(WTGAMMA-1)*TEMPERATURE",& 00489 ! error,failure) 00490 ! ELSE 00491 ! force_env%meta_env%delta_t = dummy 00492 ! ENDIF 00493 ! ELSE 00494 ! force_env%meta_env%wtgamma = 1._dp & 00495 ! + force_env%meta_env%delta_t/force_env%meta_env%wttemperature 00496 ! ENDIF 00497 ! force_env%meta_env%invdt = 1._dp/force_env%meta_env%delta_t 00498 ! ENDIF 00499 CALL tamc_force(force_env,error=error) 00500 ! CALL metadyn_write_colvar(force_env,error=error) 00501 ENDIF 00502 00503 00504 IF (simpar%do_respa)THEN 00505 CALL force_env_calc_energy_force (force_env%sub_force_env(1)%force_env,& 00506 calc_force=.TRUE.,error=error) 00507 END IF 00508 00509 ! CALL force_env_get( force_env, subsys=subsys, error=error ) 00510 ! 00511 ! CALL cp_subsys_get(subsys,atomic_kinds=atomic_kinds,local_particles=local_particles,& 00512 ! particles=particles,error=error) 00513 00514 CALL virial_evaluate(atomic_kinds%els, particles%els, local_particles,& 00515 force_env%virial, force_env%para_env%group, error=error) 00516 00517 CALL md_energy(md_env,md_ener,error) 00518 ! CALL md_write_output(md_env, error) !inits the print env at itimes == 0 also writes trajectories 00519 md_stride = 1 00520 ELSE 00521 CALL get_md_env(md_env, reftraj=reftraj, error=error) 00522 CALL initialize_reftraj(reftraj, reftraj_section, md_env, error=error) 00523 itimes = reftraj%info%first_snapshot -1 00524 md_stride = reftraj%info%stride 00525 END IF 00526 00527 CALL cp_print_key_finished_output(simpar%info_constraint, logger,& 00528 constraint_section,"CONSTRAINT_INFO",error=error) 00529 CALL cp_print_key_finished_output(simpar%lagrange_multipliers, logger,& 00530 constraint_section,"LAGRANGE_MULTIPLIERS",error=error) 00531 CALL init_mc_moves(moves) 00532 CALL init_mc_moves(gmoves) 00533 ALLOCATE (r(1:3,SIZE(particles%els)),STAT=isos) 00534 CPPostcondition(isos==0,cp_failure_level,routineP,error,failure) 00535 ! ALLOCATE (r_old(1:3,size(particles%els)),STAT=isos) 00536 ! CPPostcondition(isos==0,cp_failure_level,routineP,error,failure) 00537 CALL mc_averages_create(MCaverages) 00538 !!!!! some more buffers 00539 ! Allocate random number for Langevin Thermostat acting on COLVARS 00540 ALLOCATE (xieta(2*force_env%meta_env%n_colvar),STAT=isos) 00541 CPPostcondition(isos==0,cp_failure_level,routineP,error,failure) 00542 xieta(:) = 0.0_dp 00543 ALLOCATE (An(force_env%meta_env%n_colvar),STAT=isos) 00544 CPPostcondition(isos==0,cp_failure_level,routineP,error,failure) 00545 An(:) = 0.0_dp 00546 ALLOCATE (fz(force_env%meta_env%n_colvar),STAT=isos) 00547 CPPostcondition(isos==0,cp_failure_level,routineP,error,failure) 00548 fz(:) = 0.0_dp 00549 ALLOCATE (zbuff(2*force_env%meta_env%n_colvar),STAT=isos) 00550 CPPostcondition(isos==0,cp_failure_level,routineP,error,failure) 00551 zbuff(:) = 0.0_dp 00552 00553 IF (output_unit>0)THEN 00554 WRITE(output_unit,'(a)')"HMC|==== Initial average forces" 00555 ENDIF 00556 CALL metadyn_write_colvar_header(force_env,error=error) 00557 moves%hmc%attempts=0 00558 moves%hmc%successes=0 00559 gmoves%hmc%attempts=0 00560 gmoves%hmc%successes=0 00561 IF (initialStep==0) THEN 00562 !!! if we come from a restart we shall properly compute the average force 00563 !!! read the average force up to now 00564 DO i=1,force_env%meta_env%n_colvar 00565 fs_section => section_vals_get_subs_vals(force_env%meta_env%metadyn_section,"EXT_LAGRANGE_FS",error=error) 00566 CALL section_vals_get(fs_section, explicit=explicit, error=error) 00567 IF (explicit) THEN 00568 CALL section_vals_val_get(fs_section,"_DEFAULT_KEYWORD_",& 00569 i_rep_val=i, r_val=rval, error=error) 00570 fz(i) = rval*rand2skip 00571 END IF 00572 ENDDO 00573 00574 CALL HMCsampler(globenv,force_env,MCaverages,r,mc_par,simpar,moves,gmoves,rng_stream_mc,output_unit,& 00575 fz,zbuff,nskip=rand2skip,logger=logger,iter=0,mcsec=mc_section,mdenv=md_env,error=error) 00576 CALL cp_iterate(logger%iter_info,last=.FALSE.,iter_nr=0,error=error) 00577 CALL section_vals_val_set(mc_section,"RANDOMTOSKIP",i_val=rand2skip+nmccycles,error=error) 00578 CALL write_restart(md_env=md_env,root_section=force_env%root_section, error=error) 00579 ENDIF 00580 IF (output_unit>0)THEN 00581 WRITE(output_unit,'(a)')"HMC|==== end initial average forces" 00582 ENDIF 00583 ! call set_md_env(md_env, init=.FALSE., error=error) 00584 00585 CALL metadyn_write_colvar(force_env,error=error) 00586 00587 DO istep=1, force_env%meta_env%TAMCSteps 00588 ! Increase counters 00589 itimes = itimes + 1 00590 time = time + force_env%meta_env%dt 00591 IF (output_unit>0)THEN 00592 WRITE(output_unit,'(a)')"HMC|===================================" 00593 WRITE(output_unit,'(a,1x,i0)')"HMC| on z step ", istep 00594 ENDIF 00595 !needed when electric field fields are applied 00596 IF(ASSOCIATED(force_env%qs_env))THEN 00597 force_env%qs_env%sim_time=time 00598 force_env%qs_env%sim_step=itimes 00599 force_env%meta_env%time=force_env%qs_env%sim_time 00600 END IF 00601 00602 IF(ehrenfest_md)force_env%qs_env%rtp%istep=istep 00603 00604 CALL cp_iterate(logger%iter_info,last=(istep==force_env%meta_env%TAMCSteps),iter_nr=itimes,error=error) 00605 ! Open possible Shake output units 00606 simpar%info_constraint = cp_print_key_unit_nr(logger,constraint_section,"CONSTRAINT_INFO",& 00607 extension=".shakeLog",log_filename=.FALSE.,error=error) 00608 simpar%lagrange_multipliers = cp_print_key_unit_nr(logger,constraint_section,& 00609 "LAGRANGE_MULTIPLIERS",extension=".LagrangeMultLog",log_filename=.FALSE.,error=error) 00610 simpar%dump_lm = BTEST(cp_print_key_should_output(logger%iter_info,constraint_section,& 00611 "LAGRANGE_MULTIPLIERS",error=error),cp_p_file) 00612 00613 ! Velocity Verlet Integrator 00614 00615 moves%hmc%attempts=0 00616 moves%hmc%successes=0 00617 CALL langevinVEC(md_env,globenv,mc_env,moves,gmoves,r,& 00618 rng_stream_mc,xieta,An,fz,MCaverages,zbuff,error) 00619 00620 ! Close Shake output if requested... 00621 CALL cp_print_key_finished_output(simpar%info_constraint, logger,& 00622 constraint_section,"CONSTRAINT_INFO",error=error) 00623 CALL cp_print_key_finished_output(simpar%lagrange_multipliers, logger,& 00624 constraint_section,"LAGRANGE_MULTIPLIERS",error=error) 00625 CALL cp_iterate(logger%iter_info,iter_nr=initialStep,error=error) 00626 CALL metadyn_write_colvar(force_env,error=error) 00627 ! Free Energy calculation 00628 ! CALL free_energy_evaluate(md_env,should_stop,free_energy_section,error) 00629 00630 ![AME:UB] IF (should_stop) EXIT 00631 00632 ! Test for <PROJECT_NAME>.EXIT_MD or for WALL_TIME to exit 00633 ! Default: 00634 ! IF so we don't overwrite the restart or append to the trajectory 00635 ! because the execution could in principle stop inside the SCF where energy 00636 ! and forces are not converged. 00637 ! But: 00638 ! You can force to print the last step (for example if the method used 00639 ! to compute energy and forces is not SCF based) activating the print_key 00640 ! MOTION%MD%PRINT%FORCE_LAST. 00641 CALL external_control(should_stop,"MD",globenv=globenv,error=error) 00642 IF (should_stop) THEN 00643 CALL cp_iterate(logger%iter_info,last=.TRUE.,iter_nr=itimes,error=error) 00644 ! CALL md_output(md_env,md_section,force_env%root_section,should_stop,error=error) 00645 IF(ehrenfest_md)THEN 00646 CALL rt_prop_output(force_env%qs_env,ehrenfest,error=error) 00647 CALL rt_write_input_restart(md_env,force_env,error) 00648 END IF 00649 EXIT 00650 END IF 00651 00652 ! IF(simpar%ensemble /= reftraj_ensemble) THEN 00653 ! CALL md_energy(md_env, md_ener, error) 00654 ! CALL temperature_control(simpar, md_env, md_ener, force_env, logger, error) 00655 ! CALL comvel_control(md_ener, force_env, md_section, logger, error) 00656 ! CALL angvel_control(md_ener, force_env, md_section, logger, error) 00657 ! ELSE 00658 ! CALL md_ener_reftraj(md_env, md_ener, error) 00659 ! END IF 00660 00661 time_iter_stop=m_walltime() 00662 used_time = time_iter_stop - time_iter_start 00663 time_iter_start=time_iter_stop 00664 00665 !!!!! this writes the restart... 00666 ! CALL md_output(md_env,md_section,force_env%root_section,should_stop,error=error) 00667 00668 ! IF(simpar%ensemble == reftraj_ensemble ) THEN 00669 ! CALL write_output_reftraj(md_env,error=error) 00670 ! END IF 00671 00672 IF (output_unit>0)THEN 00673 WRITE(output_unit,'(a,1x,i0)')"HMC| end z step ", istep 00674 WRITE(output_unit,'(a)')"HMC|===================================" 00675 ENDIF 00676 END DO 00677 CALL cp_iterate(logger%iter_info,last=.TRUE.,iter_nr=itimes,error=error) 00678 force_env%qs_env%sim_time=0.0_dp 00679 force_env%qs_env%sim_step=0 00680 rand2skip=rand2skip+nmccycles*force_env%meta_env%TAMCSteps 00681 IF (initialStep == 0) rand2skip=rand2skip+nmccycles 00682 CALL section_vals_val_set(mc_section,"RANDOMTOSKIP",i_val=rand2skip,error=error) 00683 00684 CALL write_restart(md_env=md_env,root_section=force_env%root_section, error=error) 00685 ! if we need the final kinetic energy for Hybrid Monte Carlo 00686 ! hmc_ekin%final_ekin=md_ener%ekin 00687 00688 ! Remove the iteration level 00689 CALL cp_rm_iter_level(logger%iter_info,"MD",error=error) 00690 00691 ! Deallocate Thermostats and Barostats 00692 CALL release_thermostats(thermostats, error=error) 00693 CALL release_barostat_type(barostat, error=error) 00694 CALL release_simpar_type(simpar, error) 00695 CALL release_thermal_regions(thermal_regions, error) 00696 00697 CALL md_env_release(md_env, error=error) 00698 ! Clean restartable sections.. 00699 IF (my_rm_restart_info) CALL remove_restart_info(force_env%root_section,error=error) 00700 ! IF (para_env%ionode) THEN 00701 CALL delete_rng_stream(rng_stream_mc,error=error) 00702 ! ENDIF 00703 CALL MC_ENV_RELEASE(mc_env,error) 00704 DEALLOCATE(mc_par, stat=isos) ! do not forget to clean this mess 00705 CPPostcondition(isos==0,cp_failure_level,routineP,error,failure) 00706 CALL MC_MOVES_RELEASE(moves) 00707 CALL MC_MOVES_RELEASE(gmoves) 00708 DEALLOCATE(r,STAT=isos) 00709 CPPostcondition(isos==0,cp_failure_level,routineP,error,failure) 00710 ! DEALLOCATE(r_old,STAT=isos) 00711 ! CPPostcondition(isos==0,cp_failure_level,routineP,error,failure) 00712 DEALLOCATE(xieta,STAT=isos) 00713 CPPostcondition(isos==0,cp_failure_level,routineP,error,failure) 00714 DEALLOCATE(An,STAT=isos) 00715 CPPostcondition(isos==0,cp_failure_level,routineP,error,failure) 00716 DEALLOCATE(fz,STAT=isos) 00717 CPPostcondition(isos==0,cp_failure_level,routineP,error,failure) 00718 DEALLOCATE(zbuff,STAT=isos) 00719 CPPostcondition(isos==0,cp_failure_level,routineP,error,failure) 00720 CALL mc_averages_release(MCaverages) 00721 CALL timestop(handle) 00722 00723 END SUBROUTINE qs_tamc 00724 00725 00726 00727 ! ***************************************************************************** 00732 SUBROUTINE tamc_velocities_colvar(force_env,An,error) 00733 TYPE(force_env_type), POINTER :: force_env 00734 REAL(KIND=dp), DIMENSION(:), 00735 INTENT(INOUT) :: An 00736 TYPE(cp_error_type), INTENT(inout) :: error 00737 00738 CHARACTER(len=*), PARAMETER :: routineN = 'tamc_velocities_colvar', 00739 routineP = moduleN//':'//routineN 00740 00741 INTEGER :: handle, i_c 00742 LOGICAL :: failure 00743 REAL(kind=dp) :: dt, fft, sigma 00744 TYPE(cp_logger_type), POINTER :: logger 00745 TYPE(meta_env_type), POINTER :: meta_env 00746 TYPE(metavar_type), POINTER :: cv 00747 00748 failure=.FALSE. 00749 NULLIFY(logger,meta_env,cv) 00750 meta_env => force_env%meta_env 00751 CALL timeset(routineN,handle) 00752 logger => cp_error_get_logger(error) 00753 ! Add citation 00754 IF (meta_env%langevin) CALL cite_reference(VandenCic2006) 00755 dt = meta_env%dt 00756 00757 ! Evolve Velocities 00758 meta_env%epot_walls = 0.0_dp 00759 DO i_c= 1, meta_env%n_colvar 00760 cv => meta_env%metavar(i_c) 00761 fft = cv%ff_s+cv%ff_hills 00762 sigma = SQRT((meta_env%temp_wanted*kelvin)*2.0_dp*(boltzmann/joule)*cv%gamma/cv%mass) 00763 cv%vvp=cv%vvp+0.5_dp*dt*(fft/cv%mass-cv%gamma*cv%vvp)*(1.0_dp-0.25_dp*dt*cv%gamma)+An(i_c) 00764 meta_env%epot_walls = meta_env%epot_walls + cv%epot_walls 00765 ENDDO 00766 CALL timestop(handle) 00767 END SUBROUTINE tamc_velocities_colvar 00768 00769 ! ***************************************************************************** 00774 SUBROUTINE tamc_position_colvar(force_env,xieta,error) 00775 TYPE(force_env_type), POINTER :: force_env 00776 REAL(KIND=dp), DIMENSION(:), 00777 INTENT(INOUT) :: xieta 00778 TYPE(cp_error_type), INTENT(inout) :: error 00779 00780 CHARACTER(len=*), PARAMETER :: routineN = 'tamc_position_colvar', 00781 routineP = moduleN//':'//routineN 00782 00783 INTEGER :: handle, i_c 00784 LOGICAL :: failure 00785 REAL(kind=dp) :: dt, sigma 00786 TYPE(cp_logger_type), POINTER :: logger 00787 TYPE(meta_env_type), POINTER :: meta_env 00788 TYPE(metavar_type), POINTER :: cv 00789 00790 failure=.FALSE. 00791 NULLIFY(logger,meta_env,cv) 00792 meta_env => force_env%meta_env 00793 ! IF (.NOT.ASSOCIATED(meta_env)) RETURN 00794 00795 CALL timeset(routineN,handle) 00796 logger => cp_error_get_logger(error) 00797 00798 ! Add citation 00799 IF (meta_env%langevin) CALL cite_reference(VandenCic2006) 00800 dt = meta_env%dt 00801 00802 ! Update of ss0 00803 DO i_c = 1, meta_env%n_colvar 00804 cv => meta_env%metavar(i_c) 00805 sigma = SQRT((meta_env%temp_wanted*kelvin)*2.0_dp*(boltzmann/joule)*cv%gamma/cv%mass) 00806 ! cv%ss0 =cv%ss0 +dt*cv%vvp 00807 cv%ss0 =cv%ss0 +dt*cv%vvp+dt*SQRT(dt/12.0_dp)*sigma*xieta(i_c+meta_env%n_colvar) 00808 IF (cv%periodic) THEN 00809 ! A periodic COLVAR is always within [-pi,pi] 00810 cv%ss0 = SIGN(1.0_dp,ASIN(SIN(cv%ss0)))*ACOS(COS(cv%ss0)) 00811 END IF 00812 ENDDO 00813 CALL timestop(handle) 00814 00815 END SUBROUTINE tamc_position_colvar 00816 00817 ! ***************************************************************************** 00822 SUBROUTINE tamc_force(force_env,zpot,error) 00823 TYPE(force_env_type), POINTER :: force_env 00824 REAL(KIND=dp), INTENT(inout), OPTIONAL :: zpot 00825 TYPE(cp_error_type), INTENT(inout) :: error 00826 00827 CHARACTER(len=*), PARAMETER :: routineN = 'tamc_force', 00828 routineP = moduleN//':'//routineN 00829 00830 INTEGER :: handle, i, i_c, icolvar, ii 00831 LOGICAL :: explicit, failure 00832 REAL(kind=dp) :: diff_ss, dt, rval 00833 TYPE(colvar_p_type), DIMENSION(:), 00834 POINTER :: colvar_p 00835 TYPE(cp_logger_type), POINTER :: logger 00836 TYPE(cp_subsys_type), POINTER :: subsys 00837 TYPE(meta_env_type), POINTER :: meta_env 00838 TYPE(metavar_type), POINTER :: cv 00839 TYPE(particle_list_type), POINTER :: particles 00840 TYPE(section_vals_type), POINTER :: ss0_section, ss_section, 00841 vvp_section 00842 00843 failure=.FALSE. 00844 NULLIFY(logger,meta_env) 00845 meta_env => force_env%meta_env 00846 ! IF (.NOT.ASSOCIATED(meta_env)) RETURN 00847 00848 CALL timeset(routineN,handle) 00849 logger => cp_error_get_logger(error) 00850 NULLIFY(colvar_p,subsys,cv,ss0_section, vvp_section,ss_section) 00851 CALL force_env_get(force_env, subsys=subsys, error=error) 00852 00853 dt = meta_env%dt 00854 IF (.NOT.meta_env%restart) meta_env%n_steps=meta_env%n_steps+1 00855 ! compute ss and the derivative of ss with respect to the atomic positions 00856 DO i_c=1,meta_env%n_colvar 00857 cv => meta_env%metavar(i_c) 00858 icolvar = cv%icolvar 00859 CALL colvar_eval_glob_f(icolvar,force_env,error=error) 00860 cv%ss = subsys%colvar_p(icolvar)%colvar%ss 00861 ! Restart for Extended Lagrangian Metadynamics 00862 IF (meta_env%restart) THEN 00863 ! Initialize the position of the collective variable in the extended lagrange 00864 ss0_section => section_vals_get_subs_vals(meta_env%metadyn_section,"EXT_LAGRANGE_SS0",error=error) 00865 CALL section_vals_get(ss0_section, explicit=explicit, error=error) 00866 IF (explicit) THEN 00867 CALL section_vals_val_get(ss0_section,"_DEFAULT_KEYWORD_",& 00868 i_rep_val=i_c, r_val=rval, error=error) 00869 cv%ss0 = rval 00870 ELSE 00871 cv%ss0 = cv%ss 00872 END IF 00873 vvp_section => section_vals_get_subs_vals(meta_env%metadyn_section,"EXT_LAGRANGE_VVP",error=error) 00874 CALL section_vals_get(vvp_section, explicit=explicit, error=error) 00875 IF (explicit) THEN 00876 CALL setup_velocities_z(force_env,error) 00877 CALL section_vals_val_get(vvp_section,"_DEFAULT_KEYWORD_",& 00878 i_rep_val=i_c, r_val=rval, error=error) 00879 cv%vvp = rval 00880 ELSE 00881 CALL setup_velocities_z(force_env,error) 00882 ENDIF 00883 ss_section => section_vals_get_subs_vals(meta_env%metadyn_section,"EXT_LAGRANGE_SS",error=error) 00884 CALL section_vals_get(ss_section, explicit=explicit, error=error) 00885 IF (explicit) THEN 00886 CALL section_vals_val_get(ss_section,"_DEFAULT_KEYWORD_",& 00887 i_rep_val=i_c, r_val=rval, error=error) 00888 cv%ss = rval 00889 END IF 00890 END IF 00891 ! 00892 ENDDO 00893 ! forces on the atoms 00894 NULLIFY(particles) 00895 CALL cp_subsys_get(subsys, colvar_p=colvar_p, & 00896 particles=particles,error=error) 00897 00898 meta_env%restart = .FALSE. 00899 meta_env%epot_s = 0.0_dp 00900 meta_env%epot_walls = 0.0_dp 00901 DO i_c= 1, meta_env%n_colvar 00902 cv => meta_env%metavar(i_c) 00903 diff_ss = cv%ss-cv%ss0 00904 IF (cv%periodic) THEN 00905 ! The difference of a periodic COLVAR is always within [-pi,pi] 00906 diff_ss = SIGN(1.0_dp,ASIN(SIN(diff_ss)))*ACOS(COS(diff_ss)) 00907 END IF 00908 cv%epot_s = 0.5_dp*cv%lambda*diff_ss*diff_ss 00909 cv%ff_s = cv%lambda*(diff_ss) 00910 meta_env%epot_s=meta_env%epot_s+cv%epot_s 00911 icolvar = cv%icolvar 00912 00913 DO ii=1,colvar_p(icolvar)%colvar%n_atom_s 00914 i=colvar_p(icolvar)%colvar%i_atom(ii) 00915 particles%els(i)%f=particles%els(i)%f- cv%ff_s*colvar_p(icolvar)%colvar%dsdr(:,ii) 00916 ENDDO 00917 00918 ENDDO 00919 IF(PRESENT(zpot))zpot=meta_env%epot_s 00920 CALL fix_atom_control(force_env, error=error) 00921 00922 CALL timestop(handle) 00923 END SUBROUTINE tamc_force 00924 00925 00926 00927 ! ***************************************************************************** 00931 SUBROUTINE langevinVEC ( md_env, globenv,mc_env,moves,gmoves,r,& 00932 rng_stream_mc,xieta,An,fz,averages,zbuff,error) 00933 00934 TYPE(md_environment_type), POINTER :: md_env 00935 TYPE(global_environment_type), POINTER :: globenv 00936 TYPE(mc_environment_type), POINTER :: mc_env 00937 TYPE(mc_moves_type), POINTER :: moves, gmoves 00938 REAL(KIND=dp), DIMENSION(:, :), 00939 INTENT(INOUT) :: r 00940 TYPE(rng_stream_type), POINTER :: rng_stream_mc 00941 REAL(KIND=dp), DIMENSION(:), 00942 INTENT(INOUT) :: xieta, An, fz 00943 TYPE(mc_averages_type), INTENT(INOUT), 00944 POINTER :: averages 00945 REAL(KIND=dp), DIMENSION(:), 00946 INTENT(INOUT) :: zbuff 00947 TYPE(cp_error_type), INTENT(inout) :: error 00948 00949 CHARACTER(len=*), PARAMETER :: routineN = 'langevinVEC', 00950 routineP = moduleN//':'//routineN 00951 00952 INTEGER :: iprint, ivar, nparticle, 00953 nparticle_kind, nstep, 00954 output_unit 00955 LOGICAL :: failure 00956 REAL(KIND=dp) :: dt, gamma, mass, sigma 00957 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 00958 TYPE(atomic_kind_type), DIMENSION(:), 00959 POINTER :: atomic_kind_set 00960 TYPE(cell_type), POINTER :: cell 00961 TYPE(cp_logger_type), POINTER :: logger 00962 TYPE(cp_para_env_type), POINTER :: para_env 00963 TYPE(cp_subsys_type), POINTER :: subsys 00964 TYPE(distribution_1d_type), POINTER :: local_molecules, 00965 local_particles 00966 TYPE(force_env_type), POINTER :: force_env 00967 TYPE(global_constraint_type), POINTER :: gci 00968 TYPE(mc_simpar_type), POINTER :: mc_par 00969 TYPE(mol_kind_new_list_type), POINTER :: molecule_kinds 00970 TYPE(mol_new_list_type), POINTER :: molecules 00971 TYPE(molecule_kind_type), DIMENSION(:), 00972 POINTER :: molecule_kind_set 00973 TYPE(molecule_type), DIMENSION(:), 00974 POINTER :: molecule_set 00975 TYPE(particle_list_type), POINTER :: particles 00976 TYPE(particle_type), DIMENSION(:), 00977 POINTER :: particle_set 00978 TYPE(rng_stream_type), POINTER :: rng_stream 00979 TYPE(simpar_type), POINTER :: simpar 00980 TYPE(virial_type), POINTER :: virial 00981 00982 NULLIFY(logger, mc_par) 00983 logger => cp_error_get_logger(error) 00984 output_unit = cp_logger_get_default_io_unit(logger) 00985 00986 NULLIFY(rng_stream) 00987 failure = .FALSE. 00988 ! quantitites to be nullified for the get_md_env 00989 NULLIFY(simpar,force_env,para_env) 00990 ! quantities to be nullified for the force_env_get environment 00991 NULLIFY(subsys,cell) 00992 ! quantitites to be nullified for the cp_subsys_get 00993 NULLIFY(atomic_kinds,local_particles,particles,local_molecules,molecules,molecule_kinds,gci) 00994 00995 CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env,& 00996 para_env=para_env, error=error) 00997 CALL get_mc_env(mc_env, mc_par=mc_par) 00998 CALL get_mc_par(mc_par,nstep=nstep,iprint=iprint) 00999 01000 dt = simpar%dt 01001 CALL force_env_get(force_env=force_env,subsys=subsys,cell=cell,& 01002 virial=virial,error=error) 01003 01004 !!!! this bit should vanish once I understand what the hell is with it 01005 01006 ! ! Do some checks on coordinates and box 01007 CALL apply_qmmm_walls_reflective(force_env, error=error) 01008 01009 CALL cp_subsys_get(subsys=subsys,atomic_kinds=atomic_kinds,local_particles=local_particles,& 01010 particles=particles,local_molecules_new=local_molecules,molecules_new=molecules,& 01011 molecule_kinds_new=molecule_kinds,gci=gci,error=error) 01012 01013 nparticle_kind = atomic_kinds%n_els 01014 atomic_kind_set => atomic_kinds%els 01015 molecule_kind_set => molecule_kinds%els 01016 01017 nparticle = particles%n_els 01018 particle_set => particles%els 01019 molecule_set => molecules%els 01020 CPPrecondition(ASSOCIATED(force_env%meta_env),cp_failure_level,routineP,error,failure) 01021 CPPrecondition(force_env%meta_env%langevin,cp_failure_level,routineP,error,failure) 01022 ! *** Velocity Verlet for Langevin *** v(t)--> v(t+1/2) 01023 !!!!!! noise xi is in the first half, eta in the second half 01024 DO ivar = 1 , force_env%meta_env%n_colvar 01025 rng_stream => force_env%meta_env%rng(ivar)%stream 01026 xieta(ivar)=next_random_number(rng_stream,error=error) 01027 xieta(ivar+force_env%meta_env%n_colvar)=next_random_number(rng_stream,error=error) 01028 gamma=force_env%meta_env%metavar(ivar)%gamma 01029 mass=force_env%meta_env%metavar(ivar)%mass 01030 sigma = SQRT((force_env%meta_env%temp_wanted*kelvin)*2.0_dp*(boltzmann/joule)*gamma/mass) 01031 An(ivar)=0.5_dp*SQRT(dt)*sigma*(xieta(ivar)*(1.0_dp-0.5_dp*dt*gamma)-& 01032 dt*gamma*xieta(ivar+force_env%meta_env%n_colvar)/SQRT(12.0_dp)) 01033 ENDDO 01034 ! *** Velocity Verlet for Langeving *** v(t)--> v(t+1/2) 01035 CALL tamc_velocities_colvar(force_env,An,error=error) 01036 ! *** Velocity Verlet for Langevin S(t)->S(t+1) 01037 CALL tamc_position_colvar(force_env,xieta,error) 01038 !!!!! start zHMC sampler 01039 CALL HMCsampler(globenv,force_env,averages,r,mc_par,simpar,moves,gmoves,rng_stream_mc,output_unit,fz,zbuff,error=error) 01040 01041 ! CALL final_mc_write(mc_par,tmp_moves,& 01042 ! output_unit,energy_check,& 01043 ! initial_energy,final_energy,& 01044 ! averages) 01045 01046 !!!!!!!!!!!!!!!!!!!! end zHMC sampler 01047 ! *** Velocity Verlet for Langeving *** v(t+1/2)--> v(t+1) 01048 CALL tamc_velocities_colvar(force_env,An,error) 01049 ! CALL virial_evaluate ( atomic_kind_set, particle_set, & 01050 ! local_particles, virial, para_env%group, error=error) 01051 01052 END SUBROUTINE langevinVEC 01053 01054 01055 ! ***************************************************************************** 01060 01061 SUBROUTINE HMCsampler(globenv,force_env,averages,r,mc_par,simpar,moves,gmoves,rng_stream_mc,output_unit,& 01062 fz,zbuff,nskip,logger,mcsec,mdenv,iter,error) 01063 TYPE(global_environment_type), POINTER :: globenv 01064 TYPE(force_env_type), POINTER :: force_env 01065 TYPE(mc_averages_type), POINTER :: averages 01066 REAL(KIND=dp), DIMENSION(:, :), 01067 INTENT(INOUT) :: r 01068 TYPE(mc_simpar_type), POINTER :: mc_par 01069 TYPE(simpar_type), POINTER :: simpar 01070 TYPE(mc_moves_type), POINTER :: moves, gmoves 01071 TYPE(rng_stream_type), POINTER :: rng_stream_mc 01072 INTEGER, INTENT(IN) :: output_unit 01073 REAL(KIND=dp), DIMENSION(:), 01074 INTENT(INOUT) :: fz, zbuff 01075 INTEGER, INTENT(IN), OPTIONAL :: nskip 01076 TYPE(cp_logger_type), OPTIONAL, POINTER :: logger 01077 TYPE(section_vals_type), OPTIONAL, 01078 POINTER :: mcsec 01079 TYPE(md_environment_type), OPTIONAL, 01080 POINTER :: mdenv 01081 INTEGER, INTENT(IN), OPTIONAL :: iter 01082 TYPE(cp_error_type), INTENT(inout) :: error 01083 01084 INTEGER :: i, iprint, ishift, it1, j, 01085 nsamples, nstep 01086 REAL(KIND=dp) :: energy_check, old_epx, 01087 old_epz, t1 01088 TYPE(meta_env_type), POINTER :: meta_env_saved 01089 01090 IF (PRESENT(nskip)) THEN 01091 nsamples=nskip 01092 ishift=nskip 01093 ELSE 01094 nsamples=0 01095 fz=0.0_dp 01096 ishift=0 01097 ENDIF 01098 ! lrestart = .false. 01099 ! if (present(logger) .and. present(iter)) THEN 01100 ! lrestart=.true. 01101 ! ENDIF 01102 CALL get_mc_par(mc_par,nstep=nstep,iprint=iprint) 01103 meta_env_saved=> force_env%meta_env 01104 NULLIFY(force_env%meta_env) 01105 CALL force_env_get(force_env,potential_energy=old_epx,error=error) 01106 force_env%meta_env=> meta_env_saved 01107 01108 old_epz=force_env%meta_env%epot_s 01109 !!! average energy will be wrong on restarts 01110 averages%ave_energy=0.0_dp 01111 t1=force_env%qs_env%sim_time 01112 it1=force_env%qs_env%sim_step 01113 IF (output_unit > 0) THEN 01114 WRITE(output_unit,'(a,l4)')"HMC|restart? ",force_env%meta_env%restart 01115 WRITE(output_unit,'(a,3(f16.8,1x))')"HMC|Ep, Epx, Epz ",old_epx+force_env%meta_env%epot_s,old_epx,force_env%meta_env%epot_s 01116 WRITE(output_unit,'(a)')"#HMC| No | z.. | theta.. | ff_z... | ff_z/n |" 01117 ENDIF 01118 DO i = 1,nstep 01119 IF(MOD(i,iprint) == 0 .AND. (output_unit>0)) THEN 01120 WRITE(output_unit,'(a,1x,i0)') "HMC|========== On Monte Carlo cycle ",i+ishift 01121 WRITE(output_unit,'(a)') "HMC| Attempting a minitrajectory move" 01122 WRITE(output_unit,'(a,1x,i0)')"HMC| start mini-trajectory", i+ishift 01123 WRITE(output_unit,'(a,1x,i0,1x)',advance="no")"#HMC|0 ",i+ishift 01124 DO j=1,force_env%meta_env%n_colvar 01125 WRITE(output_unit,'(f16.8,1x,f16.8,1x,f16.8)',advance="no")force_env%meta_env%metavar(j)%ss0,& 01126 force_env%meta_env%metavar(j)%ss,& 01127 force_env%meta_env%metavar(j)%ff_s!,fz(j)/real(i+ishift,dp) 01128 ENDDO 01129 WRITE(output_unit,*) 01130 ENDIF 01131 01132 CALL mc_hmc_move(mc_par, force_env,globenv, moves,gmoves,old_epx,old_epz,energy_check,& 01133 r, output_unit, rng_stream_mc,zbuff,error) 01134 ! check averages... 01135 ! force average for z needed too... 01136 averages%ave_energy=averages%ave_energy*REAL(i-1,dp)/REAL(i,dp)+ 01137 old_epx/REAL(i,dp) 01138 DO j=1,force_env%meta_env%n_colvar 01139 fz(j)=fz(j)+force_env%meta_env%metavar(j)%ff_s 01140 ENDDO 01141 IF (output_unit > 0) THEN 01142 WRITE(output_unit,'(a,1x,i0)')"HMC|end mini-trajectory",i+ishift 01143 !!!!!!!! this prints z and theta(x) --ss0,ss-- needed to determine an acceptable k then the instanteneous force and some instanteneous average for force 01144 WRITE(output_unit,'(a,1x,i0,1x)',advance="no")"#HMC|1 ",i+ishift 01145 DO j=1,force_env%meta_env%n_colvar 01146 WRITE(output_unit,'(f16.8,1x,f16.8,1x,f16.8,1x,f16.8)',advance="no")force_env%meta_env%metavar(j)%ss0,& 01147 force_env%meta_env%metavar(j)%ss,& 01148 force_env%meta_env%metavar(j)%ff_s,fz(j)/REAL(i+ishift,dp) 01149 ENDDO 01150 WRITE(output_unit,*) 01151 ENDIF 01152 nsamples=nsamples+1 01153 IF(MOD(i,iprint) == 0 .AND. (output_unit>0)) THEN 01154 WRITE(output_unit,'(a,f16.8)')"HMC| Running average for potential energy ", averages%ave_energy 01155 WRITE(output_unit,'(a,1x,i0)') "HMC|======== End Monte Carlo cycle ",i+ishift 01156 ENDIF 01157 ! IF (lrestart) THEN 01158 ! k=nstep/5 01159 ! IF(MOD(i,k) == 0) THEN 01160 ! force_env%qs_env%sim_time=t1 01161 ! force_env%qs_env%sim_step=it1 01162 ! DO j=1,force_env%meta_env%n_colvar 01163 ! force_env%meta_env%metavar(j)%ff_s=fz(j)/real(i+ishift,dp) 01164 ! ENDDO 01165 ! ! CALL cp_iterate(logger%iter_info,last=.FALSE.,iter_nr=-1,error=error) 01166 ! CALL section_vals_val_set(mcsec,"RANDOMTOSKIP",i_val=i+ishift,error=error) 01167 ! CALL write_restart(md_env=mdenv,root_section=force_env%root_section, error=error) 01168 ! ! CALL cp_iterate(logger%iter_info,last=.FALSE.,iter_nr=iter,error=error) 01169 ! ENDIF 01170 ! ENDIF 01171 ENDDO 01172 force_env%qs_env%sim_time=t1 01173 force_env%qs_env%sim_step=it1 01174 IF (output_unit > 0) THEN 01175 WRITE(output_unit,'(a,i0,a,i0,a,f16.8)')"HMC| local acceptance ratio: ",moves%hmc%successes,"/" ,& 01176 moves%hmc%attempts, "=",REAL(moves%hmc%successes,dp)/ REAL(moves%hmc%attempts,dp) 01177 WRITE(output_unit,'(a,i0,a,i0,a,f16.8)')"HMC| global acceptance ratio: ",gmoves%hmc%successes,"/" ,& 01178 gmoves%hmc%attempts, "=",REAL(gmoves%hmc%successes,dp)/ REAL(gmoves%hmc%attempts,dp) 01179 ENDIF 01180 !average force 01181 DO j=1,force_env%meta_env%n_colvar 01182 force_env%meta_env%metavar(j)%ff_s=fz(j)/nsamples 01183 ENDDO 01184 END SUBROUTINE HMCsampler 01185 01186 ! ***************************************************************************** 01192 SUBROUTINE mc_hmc_move ( mc_par,force_env, globenv, moves,gmoves,old_epx,old_epz,& 01193 energy_check,r,output_unit,rng_stream,zbuff,error) 01194 01195 TYPE(mc_simpar_type), POINTER :: mc_par 01196 TYPE(force_env_type), POINTER :: force_env 01197 TYPE(global_environment_type), POINTER :: globenv 01198 TYPE(mc_moves_type), POINTER :: moves, gmoves 01199 REAL(KIND=dp), INTENT(INOUT) :: old_epx, old_epz, energy_check 01200 REAL(KIND=dp), DIMENSION(:, :), 01201 INTENT(INOUT) :: r 01202 INTEGER, INTENT(IN) :: output_unit 01203 TYPE(rng_stream_type), POINTER :: rng_stream 01204 REAL(KIND=dp), DIMENSION(:), 01205 INTENT(INOUT) :: zbuff 01206 TYPE(cp_error_type), INTENT(INOUT) :: error 01207 01208 CHARACTER(LEN=*), PARAMETER :: routineN = 'mc_hmc_move', 01209 routineP = moduleN//':'//routineN 01210 01211 INTEGER :: group, handle, iatom, istat, 01212 j, nAtoms, source 01213 LOGICAL :: failure, ionode, localise 01214 REAL(KIND=dp) :: BETA, energy_term, 01215 exp_max_val, exp_min_val, 01216 new_energy, new_epx, new_epz, 01217 rand, value, w 01218 TYPE(cp_subsys_type), POINTER :: oldsys 01219 TYPE(mc_ekin_type), POINTER :: hmc_ekin 01220 TYPE(meta_env_type), POINTER :: meta_env_saved 01221 TYPE(particle_list_type), POINTER :: particles_set 01222 TYPE(qs_scf_env_type), POINTER :: scf_env 01223 TYPE(section_vals_type), POINTER :: dft_section, input 01224 01225 ! begin the timing of the subroutine 01226 01227 CALL timeset(routineN,handle) 01228 failure=.TRUE. 01229 01230 NULLIFY(scf_env) 01231 CALL get_qs_env(force_env%qs_env,scf_env=scf_env,error=error,input=input) 01232 dft_section => section_vals_get_subs_vals(input,"DFT",error=error) 01233 01234 ! get a bunch of stuff from mc_par 01235 CALL get_mc_par(mc_par,ionode=ionode,& 01236 BETA=BETA,exp_max_val=exp_max_val,& 01237 exp_min_val=exp_min_val,source=source,group=group) 01238 01239 ! nullify some pointers 01240 ! NULLIFY(particles_set,oldsys,hmc_ekin) 01241 NULLIFY(particles_set,oldsys,meta_env_saved,hmc_ekin) 01242 ! now let's grab the particle positions 01243 CALL force_env_get(force_env,subsys=oldsys,& 01244 error=error) 01245 CALL cp_subsys_get(oldsys,particles=particles_set, & 01246 error=error) 01247 nAtoms=SIZE(particles_set%els) 01248 ! do some allocation 01249 01250 01251 ALLOCATE (hmc_ekin,STAT=istat) 01252 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01253 01254 ! record the attempt 01255 moves%hmc%attempts=moves%hmc%attempts+1 01256 gmoves%hmc%attempts=gmoves%hmc%attempts+1 01257 01258 ! save the old coordinates just in case we need to go back 01259 DO iatom=1,nAtoms 01260 r(1:3,iatom)=particles_set%els(iatom)%r(1:3) 01261 ENDDO 01262 localise=.TRUE. 01263 ! the same for collective variables data should be made,ss first half and ff_s the last half 01264 DO j=1,force_env%meta_env%n_colvar 01265 zbuff(j)=force_env%meta_env%metavar(j)%ss 01266 zbuff(j+force_env%meta_env%n_colvar)=force_env%meta_env%metavar(j)%ff_s 01267 IF ((oldsys%colvar_p(force_env%meta_env%metavar(j)%icolvar)%colvar%type_id == HBP_colvar_id) .OR.& 01268 (oldsys%colvar_p(force_env%meta_env%metavar(j)%icolvar)%colvar%type_id == WC_colvar_id)) THEN 01269 localise=.FALSE. 01270 ENDIF 01271 ENDDO 01272 01273 ! now run the MD simulation 01274 meta_env_saved=> force_env%meta_env 01275 NULLIFY(force_env%meta_env) 01276 force_env%qs_env%sim_time=0.0_dp 01277 force_env%qs_env%sim_step=0 01278 IF (.NOT. localise) THEN 01279 CALL section_vals_val_set(dft_section,"LOCALIZE%_SECTION_PARAMETERS_",l_val=.FALSE.,error=error) 01280 ENDIF 01281 CALL qs_mol_dyn(force_env,globenv,error=error,hmc_ekin=hmc_ekin) 01282 IF (.NOT. localise) THEN 01283 CALL section_vals_val_set(dft_section,"LOCALIZE%_SECTION_PARAMETERS_",l_val=.TRUE.,error=error) 01284 CALL scf_post_calculation_gpw(dft_section, scf_env, force_env%qs_env, error) 01285 ENDIF 01286 01287 CALL force_env_get(force_env, potential_energy=new_epx,error=error) 01288 01289 force_env%meta_env=>meta_env_saved 01290 CALL tamc_force(force_env, zpot=new_epz,error=error) 01291 new_energy=new_epx+new_epz 01292 IF (output_unit>0) THEN 01293 WRITE(output_unit,'(a,4(f16.8,1x))')"HMC|old Ep, Ekx, Epz, Epx ",old_epx+old_epz,hmc_ekin%initial_ekin, old_epz,old_epx 01294 WRITE(output_unit,'(a,4(f16.8,1x))')"HMC|new Ep, Ekx, Epz, Epx ",new_energy,hmc_ekin%final_ekin,new_epz,new_epx 01295 ENDIF 01296 energy_term=new_energy-old_epx-old_epz+hmc_ekin%final_ekin-hmc_ekin%initial_ekin 01297 01298 value=-BETA*(energy_term) 01299 ! to prevent overflows 01300 IF (value > exp_max_val) THEN 01301 w=10.0_dp 01302 ELSEIF(value < exp_min_val) THEN 01303 w=0.0_dp 01304 ELSE 01305 w=EXP(value) 01306 ENDIF 01307 01308 rand=next_random_number(rng_stream,error=error) 01309 IF (rand < w) THEN 01310 ! accept the move 01311 moves%hmc%successes=moves%hmc%successes+1 01312 gmoves%hmc%successes=gmoves%hmc%successes+1 01313 ! update energies 01314 energy_check=energy_check+(new_energy-old_epx-old_epz) 01315 old_epx=new_epx 01316 old_epz=new_epz 01317 ELSE 01318 ! reset the cell and particle positions 01319 DO iatom=1,nAtoms 01320 particles_set%els(iatom)%r(1:3)=r(1:3,iatom) 01321 ENDDO 01322 DO j=1,force_env%meta_env%n_colvar 01323 force_env%meta_env%metavar(j)%ss=zbuff(j) 01324 force_env%meta_env%metavar(j)%ff_s=zbuff(j+force_env%meta_env%n_colvar) 01325 ENDDO 01326 01327 ENDIF 01328 01329 DEALLOCATE(hmc_ekin,STAT=istat) 01330 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01331 01332 ! end the timing 01333 CALL timestop(handle) 01334 01335 END SUBROUTINE mc_hmc_move 01336 01337 01338 01339 SUBROUTINE metadyn_write_colvar_header(force_env,error) 01340 TYPE(force_env_type), POINTER :: force_env 01341 TYPE(cp_error_type), INTENT(inout) :: error 01342 01343 CHARACTER(len=*), PARAMETER :: routineN = 'metadyn_write_colvar_header', 01344 routineP = moduleN//':'//routineN 01345 01346 CHARACTER(len=100) :: aux, fmt 01347 CHARACTER(len=255) :: label1, label2, label3, 01348 label4, label5, label6 01349 INTEGER :: handle, i, iw, m 01350 LOGICAL :: failure 01351 TYPE(cp_logger_type), POINTER :: logger 01352 TYPE(meta_env_type), POINTER :: meta_env 01353 01354 failure=.FALSE. 01355 NULLIFY(logger,meta_env) 01356 meta_env => force_env%meta_env 01357 IF (.NOT.ASSOCIATED(meta_env)) RETURN 01358 01359 CALL timeset(routineN,handle) 01360 logger => cp_error_get_logger(error) 01361 01362 iw = cp_print_key_unit_nr(logger,meta_env%metadyn_section,& 01363 "PRINT%COLVAR",extension=".metadynLog",error=error) 01364 IF (iw>0) THEN 01365 label1="" 01366 label2="" 01367 label3="" 01368 label4="" 01369 label5="" 01370 label6="" 01371 DO i=1,meta_env%n_colvar 01372 WRITE(aux,'(a,i0)')"z_",i 01373 label1=TRIM(label1)//TRIM(aux) 01374 m=15*i-LEN_TRIM(label1)-1 01375 label1=TRIM(label1)//REPEAT(" ",m)//"|" 01376 WRITE(aux,'(a,i0)')"Theta_",i 01377 label2=TRIM(label2)//TRIM(aux) 01378 m=15*i-LEN_TRIM(label2)-1 01379 label2=TRIM(label2)//REPEAT(" ",m)//"|" 01380 WRITE(aux,'(a,i0)')"F_z",i 01381 label3=TRIM(label3)//TRIM(aux) 01382 m=15*i-LEN_TRIM(label3)-1 01383 label3=TRIM(label3)//REPEAT(" ",m)//"|" 01384 WRITE(aux,'(a,i0)')"F_h",i 01385 label4=TRIM(label4)//TRIM(aux) 01386 m=15*i-LEN_TRIM(label4)-1 01387 label4=TRIM(label4)//REPEAT(" ",m)//"|" 01388 WRITE(aux,'(a,i0)')"F_w",i 01389 label5=TRIM(label5)//TRIM(aux) 01390 m=15*i-LEN_TRIM(label5)-1 01391 label5=TRIM(label5)//REPEAT(" ",m)//"|" 01392 WRITE(aux,'(a,i0)')"v_z",i 01393 label6=TRIM(label6)//TRIM(aux) 01394 m=15*i-LEN_TRIM(label6)-1 01395 label6=TRIM(label6)//REPEAT(" ",m)//"|" 01396 ENDDO 01397 WRITE (fmt, '("(a17,6a",i0 ,",4a15)")') meta_env%n_colvar*15 01398 WRITE(iw,TRIM(fmt))"#Time[fs] |", & 01399 TRIM(label1), & 01400 TRIM(label2), & 01401 TRIM(label3), & 01402 TRIM(label4), & 01403 TRIM(label5), & 01404 TRIM(label6), & 01405 "Epot_z |", & 01406 "Ene hills |", & 01407 "Epot walls |", & 01408 "Temperature |" 01409 01410 END IF 01411 CALL cp_print_key_finished_output(iw,logger,meta_env%metadyn_section,& 01412 "PRINT%COLVAR", error=error) 01413 01414 CALL timestop(handle) 01415 01416 END SUBROUTINE metadyn_write_colvar_header 01417 01418 SUBROUTINE metadyn_write_colvar(force_env,error) 01419 TYPE(force_env_type), POINTER :: force_env 01420 TYPE(cp_error_type), INTENT(inout) :: error 01421 01422 CHARACTER(len=*), PARAMETER :: routineN = 'metadyn_write_colvar', 01423 routineP = moduleN//':'//routineN 01424 01425 INTEGER :: handle, i, i_c, iw 01426 LOGICAL :: failure 01427 REAL(KIND=dp) :: temp 01428 TYPE(cp_logger_type), POINTER :: logger 01429 TYPE(meta_env_type), POINTER :: meta_env 01430 TYPE(metavar_type), POINTER :: cv 01431 01432 failure=.FALSE. 01433 NULLIFY(logger,meta_env,cv) 01434 meta_env => force_env%meta_env 01435 IF (.NOT.ASSOCIATED(meta_env)) RETURN 01436 01437 CALL timeset(routineN,handle) 01438 logger => cp_error_get_logger(error) 01439 01440 IF (meta_env%langevin) THEN 01441 meta_env%ekin_s = 0.0_dp 01442 ! meta_env%epot_s = 0.0_dp 01443 DO i_c= 1, meta_env%n_colvar 01444 cv => meta_env%metavar(i_c) 01445 meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2 01446 ENDDO 01447 END IF 01448 01449 ! write COLVAR file 01450 iw = cp_print_key_unit_nr(logger,meta_env%metadyn_section,& 01451 "PRINT%COLVAR",extension=".metadynLog",error=error) 01452 IF (iw>0) THEN 01453 IF (meta_env%extended_lagrange) THEN 01454 WRITE(iw,'(f16.8,70f15.8)')meta_env%time*femtoseconds, & 01455 (meta_env%metavar(i)%ss0,i=1,meta_env%n_colvar), & 01456 (meta_env%metavar(i)%ss,i=1,meta_env%n_colvar), & 01457 (meta_env%metavar(i)%ff_s,i=1,meta_env%n_colvar), & 01458 (meta_env%metavar(i)%ff_hills,i=1,meta_env%n_colvar), & 01459 (meta_env%metavar(i)%ff_walls,i=1,meta_env%n_colvar), & 01460 (meta_env%metavar(i)%vvp,i=1,meta_env%n_colvar), & 01461 meta_env%epot_s, & 01462 meta_env%hills_env%energy, & 01463 meta_env%epot_walls, & 01464 (meta_env%ekin_s)*2.0_dp/(REAL(meta_env%n_colvar,KIND=dp))*kelvin 01465 ELSE 01466 WRITE(iw,'(f16.8,40f13.5)') meta_env%time*femtoseconds,& 01467 (meta_env%metavar(i)%ss0,i=1,meta_env%n_colvar),& 01468 (meta_env%metavar(i)%ff_hills,i=1,meta_env%n_colvar),& 01469 (meta_env%metavar(i)%ff_walls,i=1,meta_env%n_colvar),& 01470 meta_env%hills_env%energy,& 01471 meta_env%epot_walls 01472 END IF 01473 END IF 01474 CALL cp_print_key_finished_output(iw,logger,meta_env%metadyn_section,& 01475 "PRINT%COLVAR", error=error) 01476 ! Temperature for COLVAR 01477 IF (meta_env%extended_lagrange) THEN 01478 temp = meta_env%ekin_s*2.0_dp/(REAL(meta_env%n_colvar,KIND=dp))*kelvin 01479 meta_env%avg_temp = (meta_env%avg_temp*REAL(meta_env%n_steps,KIND=dp)+ 01480 temp)/REAL(meta_env%n_steps+1,KIND=dp) 01481 iw = cp_print_key_unit_nr(logger,meta_env%metadyn_section,& 01482 "PRINT%TEMPERATURE_COLVAR",extension=".metadynLog",error=error) 01483 IF (iw > 0) THEN 01484 WRITE (iw, '(T2,79("-"))') 01485 WRITE (iw,'( A,T51,f10.2,T71,f10.2)' )' COLVARS INSTANTANEOUS/AVERAGE TEMPERATURE ',& 01486 temp, meta_env%avg_temp 01487 WRITE (iw, '(T2,79("-"))') 01488 ENDIF 01489 CALL cp_print_key_finished_output(iw,logger,meta_env%metadyn_section,& 01490 "PRINT%TEMPERATURE_COLVAR", error=error) 01491 END IF 01492 CALL timestop(handle) 01493 01494 END SUBROUTINE metadyn_write_colvar 01495 01496 01497 SUBROUTINE setup_velocities_z(force_env,error) 01498 TYPE(force_env_type), POINTER :: force_env 01499 TYPE(cp_error_type), INTENT(inout) :: error 01500 01501 INTEGER :: i_c 01502 REAL(kind=dp) :: ekin_w, fac_t 01503 TYPE(meta_env_type), POINTER :: meta_env 01504 TYPE(metavar_type), POINTER :: cv 01505 01506 NULLIFY(meta_env) 01507 meta_env=>force_env%meta_env 01508 meta_env%ekin_s = 0.0_dp 01509 DO i_c=1,meta_env%n_colvar 01510 cv => meta_env%metavar(i_c) 01511 cv%vvp = next_random_number(force_env%globenv%gaussian_rng_stream,error=error) 01512 meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2 01513 END DO 01514 ekin_w = 0.5_dp*meta_env%temp_wanted*REAL(meta_env%n_colvar,KIND=dp) 01515 fac_t = SQRT(ekin_w/MAX(meta_env%ekin_s,1.0E-8_dp)) 01516 DO i_c = 1, meta_env%n_colvar 01517 cv => meta_env%metavar(i_c) 01518 cv%vvp = cv%vvp*fac_t 01519 ENDDO 01520 END SUBROUTINE setup_velocities_z 01521 END MODULE tamc_run
1.7.3