|
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 ! ***************************************************************************** 00013 MODULE extended_system_init 00014 00015 USE cell_types, ONLY: cell_type 00016 USE cp_para_types, ONLY: cp_para_env_type 00017 USE distribution_1d_types, ONLY: distribution_1d_type 00018 USE extended_system_mapping, ONLY: nhc_to_barostat_mapping,& 00019 nhc_to_particle_mapping,& 00020 nhc_to_particle_mapping_fast,& 00021 nhc_to_particle_mapping_slow,& 00022 nhc_to_shell_mapping 00023 USE extended_system_types, ONLY: debug_isotropic_limit,& 00024 lnhc_parameters_type,& 00025 map_info_type,& 00026 npt_info_type 00027 USE f77_blas 00028 USE global_types, ONLY: global_environment_type 00029 USE input_constants, ONLY: do_thermo_only_master,& 00030 npe_f_ensemble,& 00031 npe_i_ensemble,& 00032 nph_uniaxial_damped_ensemble,& 00033 nph_uniaxial_ensemble,& 00034 npt_f_ensemble,& 00035 npt_i_ensemble 00036 USE input_cp2k_binary_restarts, ONLY: read_binary_thermostats_nose 00037 USE input_section_types, ONLY: section_vals_get,& 00038 section_vals_get_subs_vals,& 00039 section_vals_remove_values,& 00040 section_vals_type,& 00041 section_vals_val_get 00042 USE kinds, ONLY: dp 00043 USE molecule_kind_types, ONLY: molecule_kind_type 00044 USE molecule_types_new, ONLY: global_constraint_type,& 00045 molecule_type 00046 USE parallel_rng_types, ONLY: next_random_number 00047 USE simpar_types, ONLY: simpar_type 00048 USE termination, ONLY: stop_program 00049 USE thermostat_types, ONLY: thermostat_info_type 00050 USE thermostat_utils, ONLY: get_nhc_energies 00051 USE timings, ONLY: timeset,& 00052 timestop 00053 #include "cp_common_uses.h" 00054 00055 IMPLICIT NONE 00056 00057 PRIVATE 00058 00059 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'extended_system_init' 00060 00061 PUBLIC :: initialize_nhc_part, initialize_nhc_baro, initialize_npt, & 00062 initialize_nhc_shell, initialize_nhc_slow, initialize_nhc_fast 00063 00064 CONTAINS 00065 00066 ! ***************************************************************************** 00069 SUBROUTINE initialize_npt ( simpar, globenv, npt_info, cell, work_section, error) 00070 00071 TYPE(simpar_type), POINTER :: simpar 00072 TYPE(global_environment_type), POINTER :: globenv 00073 TYPE(npt_info_type), DIMENSION(:, :), 00074 POINTER :: npt_info 00075 TYPE(cell_type), POINTER :: cell 00076 TYPE(section_vals_type), POINTER :: work_section 00077 TYPE(cp_error_type), INTENT(inout) :: error 00078 00079 CHARACTER(LEN=*), PARAMETER :: routineN = 'initialize_npt', 00080 routineP = moduleN//':'//routineN 00081 00082 INTEGER :: handle, i, ind, j, stat 00083 LOGICAL :: explicit, failure, restart 00084 REAL(KIND=dp) :: temp 00085 REAL(KIND=dp), DIMENSION(:), POINTER :: buffer 00086 TYPE(section_vals_type), POINTER :: work_section2 00087 00088 CALL timeset(routineN,handle) 00089 00090 NULLIFY (work_section2) 00091 00092 explicit = .FALSE. 00093 failure=.FALSE. 00094 restart = .FALSE. 00095 00096 CPPrecondition(.NOT.ASSOCIATED(npt_info),cp_fatal_level,routineP,error,failure) 00097 00098 ! first allocating the npt_info_type if requested 00099 SELECT CASE ( simpar % ensemble ) 00100 CASE ( npt_i_ensemble, npe_i_ensemble ) 00101 ALLOCATE ( npt_info ( 1, 1 ), STAT = stat ) 00102 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00103 npt_info ( :, : ) % eps = LOG ( cell % deth ) / 3.0_dp 00104 temp = simpar % temp_baro_ext 00105 00106 CASE ( npt_f_ensemble, npe_f_ensemble ) 00107 ALLOCATE ( npt_info ( 3, 3 ), STAT = stat ) 00108 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00109 temp = simpar % temp_baro_ext 00110 00111 CASE ( nph_uniaxial_ensemble ) 00112 ALLOCATE ( npt_info ( 1, 1 ), STAT = stat ) 00113 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00114 temp = simpar% temp_baro_ext 00115 00116 CASE ( nph_uniaxial_damped_ensemble ) 00117 ALLOCATE ( npt_info ( 1, 1 ), STAT = stat ) 00118 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00119 temp = simpar % temp_baro_ext 00120 00121 CASE DEFAULT 00122 ! Do nothing.. 00123 NULLIFY(npt_info) 00124 END SELECT 00125 00126 IF (ASSOCIATED(npt_info)) THEN 00127 IF (ASSOCIATED(work_section)) THEN 00128 work_section2 => section_vals_get_subs_vals(work_section,"VELOCITY",error=error) 00129 CALL section_vals_get(work_section2, explicit=explicit, error=error) 00130 restart=explicit 00131 work_section2 => section_vals_get_subs_vals(work_section,"MASS",error=error) 00132 CALL section_vals_get(work_section2, explicit=explicit, error=error) 00133 CALL cp_assert(restart.EQV.explicit,cp_failure_level,cp_assertion_failed,routineP,& 00134 "You need to define both VELOCITY and MASS section (or none) in the BAROSTAT section",& 00135 error=error, failure=failure) 00136 restart=explicit.AND.restart 00137 END IF 00138 00139 IF ( restart ) THEN 00140 CALL section_vals_val_get(work_section,"VELOCITY%_DEFAULT_KEYWORD_",r_vals=buffer,error=error) 00141 ind = 0 00142 DO i = 1, SIZE(npt_info,1) 00143 DO j = 1, SIZE(npt_info,2) 00144 ind = ind + 1 00145 npt_info ( i, j ) % v = buffer(ind) 00146 END DO 00147 END DO 00148 CALL section_vals_val_get(work_section,"MASS%_DEFAULT_KEYWORD_",r_vals=buffer,error=error) 00149 ind = 0 00150 DO i = 1, SIZE(npt_info,1) 00151 DO j = 1, SIZE(npt_info,2) 00152 ind = ind + 1 00153 npt_info ( i, j ) % mass = buffer(ind) 00154 END DO 00155 END DO 00156 ELSE 00157 CALL init_barostat_variables ( npt_info, simpar % tau_cell, temp, & 00158 simpar % nfree, simpar % ensemble, simpar % cmass,& 00159 globenv ,error=error) 00160 END IF 00161 00162 END IF 00163 00164 CALL timestop(handle) 00165 00166 END SUBROUTINE initialize_npt 00167 00168 ! ***************************************************************************** 00172 SUBROUTINE initialize_nhc_baro ( simpar, para_env, globenv, nhc, nose_section, save_mem, error ) 00173 00174 TYPE(simpar_type), POINTER :: simpar 00175 TYPE(cp_para_env_type), POINTER :: para_env 00176 TYPE(global_environment_type), POINTER :: globenv 00177 TYPE(lnhc_parameters_type), POINTER :: nhc 00178 TYPE(section_vals_type), POINTER :: nose_section 00179 LOGICAL, INTENT(IN) :: save_mem 00180 TYPE(cp_error_type), INTENT(inout) :: error 00181 00182 CHARACTER(len=*), PARAMETER :: routineN = 'initialize_nhc_baro', 00183 routineP = moduleN//':'//routineN 00184 00185 INTEGER :: handle, stat 00186 LOGICAL :: failure, restart 00187 REAL(KIND=dp) :: temp 00188 00189 CALL timeset(routineN,handle) 00190 00191 failure = .FALSE. 00192 restart = .FALSE. 00193 00194 CALL nhc_to_barostat_mapping ( simpar, nhc, error) 00195 00196 ! Set up the Yoshida weights 00197 IF ( nhc % nyosh > 0 ) THEN 00198 ALLOCATE ( nhc % dt_yosh ( 1 : nhc%nyosh ), STAT = stat ) 00199 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00200 CALL set_yoshida_coef ( nhc, simpar % dt ) 00201 END IF 00202 00203 CALL restart_nose(nhc,nose_section,save_mem,restart,"","",para_env,error) 00204 00205 IF (.NOT.restart) THEN 00206 ! Initializing thermostat forces and velocities for the Nose-Hoover 00207 ! Chain variables 00208 SELECT CASE (simpar % ensemble ) 00209 CASE DEFAULT 00210 temp = simpar%temp_baro_ext 00211 END SELECT 00212 IF ( nhc % nhc_len /= 0 ) THEN 00213 CALL init_nhc_variables(nhc,temp,para_env,globenv,error=error) 00214 END IF 00215 END IF 00216 00217 CALL init_nhc_forces ( nhc, error=error) 00218 00219 CALL timestop(handle) 00220 00221 END SUBROUTINE initialize_nhc_baro 00222 00223 ! ***************************************************************************** 00226 SUBROUTINE initialize_nhc_slow ( thermostat_info, simpar, local_molecules,& 00227 molecule, molecule_kind_set, para_env, globenv, nhc, nose_section,& 00228 gci, save_mem, error) 00229 00230 TYPE(thermostat_info_type), POINTER :: thermostat_info 00231 TYPE(simpar_type), POINTER :: simpar 00232 TYPE(distribution_1d_type), POINTER :: local_molecules 00233 TYPE(molecule_type), POINTER :: molecule( : ) 00234 TYPE(molecule_kind_type), POINTER :: molecule_kind_set( : ) 00235 TYPE(cp_para_env_type), POINTER :: para_env 00236 TYPE(global_environment_type), POINTER :: globenv 00237 TYPE(lnhc_parameters_type), POINTER :: nhc 00238 TYPE(section_vals_type), POINTER :: nose_section 00239 TYPE(global_constraint_type), POINTER :: gci 00240 LOGICAL, INTENT(IN) :: save_mem 00241 TYPE(cp_error_type), INTENT(INOUT) :: error 00242 00243 CHARACTER(len=*), PARAMETER :: routineN = 'initialize_nhc_slow', 00244 routineP = moduleN//':'//routineN 00245 00246 INTEGER :: handle, stat 00247 LOGICAL :: failure, restart 00248 00249 CALL timeset(routineN,handle) 00250 00251 failure = .FALSE. 00252 restart = .FALSE. 00253 ! fire up the thermostats, if not NVE 00254 00255 CALL nhc_to_particle_mapping_slow ( thermostat_info, simpar, local_molecules,& 00256 molecule, molecule_kind_set, nhc, para_env, gci, error ) 00257 00258 ! Set up the Yoshida weights 00259 IF ( nhc % nyosh > 0 ) THEN 00260 ALLOCATE ( nhc % dt_yosh ( 1 : nhc%nyosh ), STAT = stat ) 00261 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00262 CALL set_yoshida_coef ( nhc, simpar % dt) 00263 END IF 00264 00265 CALL restart_nose(nhc,nose_section,save_mem,restart,"","",para_env,error) 00266 00267 IF (.NOT.restart) THEN 00268 ! Initializing thermostat forces and velocities for the Nose-Hoover 00269 ! Chain variables 00270 IF ( nhc % nhc_len /= 0 ) THEN 00271 CALL init_nhc_variables(nhc,simpar%temp_slow,para_env,globenv,error=error) 00272 END IF 00273 END IF 00274 00275 CALL init_nhc_forces ( nhc, error=error) 00276 00277 CALL timestop(handle) 00278 00279 END SUBROUTINE initialize_nhc_slow 00280 00281 ! ***************************************************************************** 00284 SUBROUTINE initialize_nhc_fast ( thermostat_info, simpar, local_molecules,& 00285 molecule, molecule_kind_set, para_env, globenv, nhc, nose_section,& 00286 gci, save_mem, error) 00287 00288 TYPE(thermostat_info_type), POINTER :: thermostat_info 00289 TYPE(simpar_type), POINTER :: simpar 00290 TYPE(distribution_1d_type), POINTER :: local_molecules 00291 TYPE(molecule_type), POINTER :: molecule( : ) 00292 TYPE(molecule_kind_type), POINTER :: molecule_kind_set( : ) 00293 TYPE(cp_para_env_type), POINTER :: para_env 00294 TYPE(global_environment_type), POINTER :: globenv 00295 TYPE(lnhc_parameters_type), POINTER :: nhc 00296 TYPE(section_vals_type), POINTER :: nose_section 00297 TYPE(global_constraint_type), POINTER :: gci 00298 LOGICAL, INTENT(IN) :: save_mem 00299 TYPE(cp_error_type), INTENT(INOUT) :: error 00300 00301 CHARACTER(len=*), PARAMETER :: routineN = 'initialize_nhc_fast', 00302 routineP = moduleN//':'//routineN 00303 00304 INTEGER :: handle, stat 00305 LOGICAL :: failure, restart 00306 00307 CALL timeset(routineN,handle) 00308 00309 failure = .FALSE. 00310 restart = .FALSE. 00311 ! fire up the thermostats, if not NVE 00312 00313 CALL nhc_to_particle_mapping_fast ( thermostat_info, simpar, local_molecules,& 00314 molecule, molecule_kind_set, nhc, para_env, gci, error ) 00315 00316 ! Set up the Yoshida weights 00317 IF ( nhc % nyosh > 0 ) THEN 00318 ALLOCATE ( nhc % dt_yosh ( 1 : nhc%nyosh ), STAT = stat ) 00319 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00320 CALL set_yoshida_coef ( nhc, simpar % dt) 00321 END IF 00322 00323 CALL restart_nose(nhc,nose_section,save_mem,restart,"","",para_env,error) 00324 00325 IF (.NOT.restart) THEN 00326 ! Initializing thermostat forces and velocities for the Nose-Hoover 00327 ! Chain variables 00328 IF ( nhc % nhc_len /= 0 ) THEN 00329 CALL init_nhc_variables(nhc,simpar%temp_fast,para_env,globenv,error=error) 00330 END IF 00331 END IF 00332 00333 CALL init_nhc_forces ( nhc, error=error) 00334 00335 CALL timestop(handle) 00336 00337 END SUBROUTINE initialize_nhc_fast 00338 00339 ! ***************************************************************************** 00342 SUBROUTINE initialize_nhc_part ( thermostat_info, simpar, local_molecules,& 00343 molecule, molecule_kind_set, para_env, globenv, nhc, nose_section,& 00344 gci, save_mem, binary_restart_file_name, error) 00345 00346 TYPE(thermostat_info_type), POINTER :: thermostat_info 00347 TYPE(simpar_type), POINTER :: simpar 00348 TYPE(distribution_1d_type), POINTER :: local_molecules 00349 TYPE(molecule_type), POINTER :: molecule( : ) 00350 TYPE(molecule_kind_type), POINTER :: molecule_kind_set( : ) 00351 TYPE(cp_para_env_type), POINTER :: para_env 00352 TYPE(global_environment_type), POINTER :: globenv 00353 TYPE(lnhc_parameters_type), POINTER :: nhc 00354 TYPE(section_vals_type), POINTER :: nose_section 00355 TYPE(global_constraint_type), POINTER :: gci 00356 LOGICAL, INTENT(IN) :: save_mem 00357 CHARACTER(LEN=*), INTENT(IN) :: binary_restart_file_name 00358 TYPE(cp_error_type), INTENT(INOUT) :: error 00359 00360 CHARACTER(len=*), PARAMETER :: routineN = 'initialize_nhc_part', 00361 routineP = moduleN//':'//routineN 00362 00363 INTEGER :: handle, stat 00364 LOGICAL :: failure, restart 00365 00366 CALL timeset(routineN,handle) 00367 00368 failure = .FALSE. 00369 restart = .FALSE. 00370 ! fire up the thermostats, if not NVE 00371 00372 CALL nhc_to_particle_mapping ( thermostat_info, simpar, local_molecules,& 00373 molecule, molecule_kind_set, nhc, para_env, gci, error ) 00374 00375 ! Set up the Yoshida weights 00376 IF ( nhc % nyosh > 0 ) THEN 00377 ALLOCATE ( nhc % dt_yosh ( 1 : nhc%nyosh ), STAT = stat ) 00378 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00379 CALL set_yoshida_coef ( nhc, simpar % dt) 00380 END IF 00381 00382 CALL restart_nose(nhc,nose_section,save_mem,restart,binary_restart_file_name,& 00383 "PARTICLE",para_env,error) 00384 00385 IF (.NOT.restart) THEN 00386 ! Initializing thermostat forces and velocities for the Nose-Hoover 00387 ! Chain variables 00388 IF ( nhc % nhc_len /= 0 ) THEN 00389 CALL init_nhc_variables(nhc,simpar%temp_ext,para_env,globenv,error=error) 00390 END IF 00391 END IF 00392 00393 CALL init_nhc_forces ( nhc, error=error) 00394 00395 CALL timestop(handle) 00396 00397 END SUBROUTINE initialize_nhc_part 00398 00399 ! ***************************************************************************** 00402 SUBROUTINE initialize_nhc_shell( thermostat_info, simpar, local_molecules,& 00403 molecule, molecule_kind_set, para_env, globenv, nhc, nose_section,& 00404 gci, save_mem, binary_restart_file_name, error) 00405 00406 TYPE(thermostat_info_type), POINTER :: thermostat_info 00407 TYPE(simpar_type), POINTER :: simpar 00408 TYPE(distribution_1d_type), POINTER :: local_molecules 00409 TYPE(molecule_type), POINTER :: molecule( : ) 00410 TYPE(molecule_kind_type), POINTER :: molecule_kind_set( : ) 00411 TYPE(cp_para_env_type), POINTER :: para_env 00412 TYPE(global_environment_type), POINTER :: globenv 00413 TYPE(lnhc_parameters_type), POINTER :: nhc 00414 TYPE(section_vals_type), POINTER :: nose_section 00415 TYPE(global_constraint_type), POINTER :: gci 00416 LOGICAL, INTENT(IN) :: save_mem 00417 CHARACTER(LEN=*), INTENT(IN) :: binary_restart_file_name 00418 TYPE(cp_error_type), INTENT(inout) :: error 00419 00420 CHARACTER(len=*), PARAMETER :: routineN = 'initialize_nhc_shell', 00421 routineP = moduleN//':'//routineN 00422 00423 INTEGER :: handle, stat 00424 LOGICAL :: failure, restart 00425 00426 CALL timeset(routineN,handle) 00427 00428 failure = .FALSE. 00429 CALL nhc_to_shell_mapping(thermostat_info, simpar, local_molecules,& 00430 molecule, molecule_kind_set, nhc, para_env, gci, error) 00431 00432 restart = .FALSE. 00433 ! Set up the Yoshida weights 00434 IF ( nhc % nyosh > 0 ) THEN 00435 ALLOCATE ( nhc % dt_yosh ( 1 : nhc%nyosh ), STAT = stat ) 00436 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00437 CALL set_yoshida_coef ( nhc, simpar % dt ) 00438 END IF 00439 00440 CALL restart_nose(nhc,nose_section,save_mem,restart,binary_restart_file_name,& 00441 "SHELL",para_env,error) 00442 00443 IF (.NOT.restart) THEN 00444 ! Initialize thermostat forces and velocities 00445 ! Chain variables 00446 IF ( nhc % nhc_len /= 0 ) THEN 00447 CALL init_nhc_variables (nhc, simpar%temp_sh_ext, para_env, globenv, error=error) 00448 END IF 00449 END IF 00450 00451 CALL init_nhc_forces ( nhc, error=error) 00452 00453 CALL timestop(handle) 00454 00455 END SUBROUTINE initialize_nhc_shell 00456 00457 ! ***************************************************************************** 00464 SUBROUTINE set_yoshida_coef ( nhc, dt ) 00465 00466 TYPE(lnhc_parameters_type), POINTER :: nhc 00467 REAL(KIND=dp), INTENT(IN) :: dt 00468 00469 CHARACTER(len=*), PARAMETER :: routineN = 'set_yoshida_coef', 00470 routineP = moduleN//':'//routineN 00471 00472 REAL(KIND=dp), DIMENSION(nhc%nyosh) :: yosh_wt 00473 00474 SELECT CASE (nhc % nyosh) 00475 CASE DEFAULT 00476 CALL stop_program(routineN,moduleN,__LINE__,'Value not available.') 00477 CASE (1) 00478 yosh_wt(1) = 1.0_dp 00479 CASE (3) 00480 yosh_wt(1) = 1.0_dp/(2.0_dp-(2.0_dp)**(1.0_dp/3.0_dp)) 00481 yosh_wt(2) = 1.0_dp - 2.0_dp*yosh_wt(1) 00482 yosh_wt(3) = yosh_wt(1) 00483 CASE (5) 00484 yosh_wt(1) = 1.0_dp/(4.0_dp-(4.0_dp)**(1.0_dp/3.0_dp)) 00485 yosh_wt(2) = yosh_wt(1) 00486 yosh_wt(4) = yosh_wt(1) 00487 yosh_wt(5) = yosh_wt(1) 00488 yosh_wt(3) = 1.0_dp - 4.0_dp*yosh_wt(1) 00489 CASE (7) 00490 yosh_wt(1) = .78451361047756_dp 00491 yosh_wt(2) = .235573213359357_dp 00492 yosh_wt(3) = -1.17767998417887_dp 00493 yosh_wt(4) = 1.0_dp - 2.0_dp*(yosh_wt(1)+yosh_wt(2)+yosh_wt(3)) 00494 yosh_wt(5) = yosh_wt(3) 00495 yosh_wt(6) = yosh_wt(2) 00496 yosh_wt(7) = yosh_wt(1) 00497 CASE (9) 00498 yosh_wt(1) = 0.192_dp 00499 yosh_wt(2) = 0.554910818409783619692725006662999_dp 00500 yosh_wt(3) = 0.124659619941888644216504240951585_dp 00501 yosh_wt(4) = -0.843182063596933505315033808282941_dp 00502 yosh_wt(5) = 1.0_dp - 2.0_dp*(yosh_wt(1)+yosh_wt(2)+& 00503 yosh_wt(3)+yosh_wt(4)) 00504 yosh_wt(6) = yosh_wt(4) 00505 yosh_wt(7) = yosh_wt(3) 00506 yosh_wt(8) = yosh_wt(2) 00507 yosh_wt(9) = yosh_wt(1) 00508 CASE (15) 00509 yosh_wt(1) = 0.102799849391985_dp 00510 yosh_wt(2) = -0.196061023297549e1_dp 00511 yosh_wt(3) = 0.193813913762276e1_dp 00512 yosh_wt(4) = -0.158240635368243_dp 00513 yosh_wt(5) = -0.144485223686048e1_dp 00514 yosh_wt(6) = 0.253693336566229_dp 00515 yosh_wt(7) = 0.914844246229740_dp 00516 yosh_wt(8) = 1.0_dp - 2.0_dp*(yosh_wt(1)+yosh_wt(2)+& 00517 yosh_wt(3)+yosh_wt(4)+yosh_wt(5)+yosh_wt(6)+yosh_wt(7)) 00518 yosh_wt(9) = yosh_wt(7) 00519 yosh_wt(10) = yosh_wt(6) 00520 yosh_wt(11) = yosh_wt(5) 00521 yosh_wt(12) = yosh_wt(4) 00522 yosh_wt(13) = yosh_wt(3) 00523 yosh_wt(14) = yosh_wt(2) 00524 yosh_wt(15) = yosh_wt(1) 00525 END SELECT 00526 nhc % dt_yosh = dt * yosh_wt / REAL ( nhc % nc,KIND=dp) 00527 00528 END SUBROUTINE set_yoshida_coef 00529 00530 ! ***************************************************************************** 00537 SUBROUTINE restart_nose(nhc,nose_section,save_mem,restart,& 00538 binary_restart_file_name,thermostat_name,& 00539 para_env,error) 00540 00541 TYPE(lnhc_parameters_type), POINTER :: nhc 00542 TYPE(section_vals_type), POINTER :: nose_section 00543 LOGICAL, INTENT(IN) :: save_mem 00544 LOGICAL, INTENT(OUT) :: restart 00545 CHARACTER(LEN=*), INTENT(IN) :: binary_restart_file_name, 00546 thermostat_name 00547 TYPE(cp_para_env_type), POINTER :: para_env 00548 TYPE(cp_error_type), INTENT(inout) :: error 00549 00550 CHARACTER(len=*), PARAMETER :: routineN = 'restart_nose', 00551 routineP = moduleN//':'//routineN 00552 00553 INTEGER :: handle, i, ind, j 00554 LOGICAL :: explicit 00555 REAL(KIND=dp), DIMENSION(:), POINTER :: buffer 00556 TYPE(map_info_type), POINTER :: map_info 00557 TYPE(section_vals_type), POINTER :: work_section 00558 00559 CALL timeset(routineN,handle) 00560 00561 NULLIFY (buffer) 00562 NULLIFY (work_section) 00563 00564 IF (LEN_TRIM(binary_restart_file_name) > 0) THEN 00565 00566 ! Read binary restart file, if available 00567 00568 CALL read_binary_thermostats_nose(thermostat_name,nhc,binary_restart_file_name,& 00569 restart,para_env,error) 00570 00571 ELSE 00572 00573 ! Read the default restart file in ASCII format 00574 00575 explicit = .FALSE. 00576 restart = .FALSE. 00577 00578 IF (ASSOCIATED(nose_section)) THEN 00579 work_section => section_vals_get_subs_vals(nose_section,"VELOCITY",error=error) 00580 CALL section_vals_get(work_section, explicit=explicit, error=error) 00581 restart=explicit 00582 work_section => section_vals_get_subs_vals(nose_section,"COORD",error=error) 00583 CALL section_vals_get(work_section, explicit=explicit, error=error) 00584 CALL cp_assert(restart.or..not.explicit,cp_failure_level,cp_assertion_failed,routineP,& 00585 "You need to define both VELOCITY and COORD and MASS and FORCE section (or none) in the NOSE section",& 00586 error=error) 00587 restart=explicit.and.restart 00588 work_section => section_vals_get_subs_vals(nose_section,"MASS",error=error) 00589 CALL section_vals_get(work_section, explicit=explicit, error=error) 00590 CALL cp_assert(restart.or..not.explicit,cp_failure_level,cp_assertion_failed,routineP,& 00591 "You need to define both VELOCITY and COORD and MASS and FORCE section (or none) in the NOSE section",& 00592 error=error) 00593 restart=explicit.and.restart 00594 work_section => section_vals_get_subs_vals(nose_section,"FORCE",error=error) 00595 CALL section_vals_get(work_section, explicit=explicit, error=error) 00596 CALL cp_assert(restart.or..not.explicit,cp_failure_level,cp_assertion_failed,routineP,& 00597 "You need to define both VELOCITY and COORD and MASS and FORCE section (or none) in the NOSE section",& 00598 error=error) 00599 restart=explicit.and.restart 00600 END IF 00601 00602 IF (restart) THEN 00603 map_info => nhc%map_info 00604 CALL section_vals_val_get(nose_section,"COORD%_DEFAULT_KEYWORD_",r_vals=buffer,error=error) 00605 DO i = 1, SIZE ( nhc % nvt, 2) 00606 ind = map_info%index(i) 00607 ind = (ind-1) * nhc % nhc_len 00608 DO j = 1, SIZE(nhc % nvt, 1) 00609 ind = ind + 1 00610 nhc % nvt(j,i) % eta = buffer ( ind ) 00611 END DO 00612 END DO 00613 CALL section_vals_val_get(nose_section,"VELOCITY%_DEFAULT_KEYWORD_",r_vals=buffer,error=error) 00614 DO i = 1, SIZE ( nhc % nvt, 2) 00615 ind = map_info%index(i) 00616 ind = (ind-1) * nhc % nhc_len 00617 DO j = 1, SIZE(nhc % nvt, 1) 00618 ind = ind + 1 00619 nhc % nvt(j,i) % v = buffer ( ind ) 00620 END DO 00621 END DO 00622 CALL section_vals_val_get(nose_section,"MASS%_DEFAULT_KEYWORD_",r_vals=buffer,error=error) 00623 DO i = 1, SIZE ( nhc % nvt, 2) 00624 ind = map_info%index(i) 00625 ind = (ind-1) * nhc % nhc_len 00626 DO j = 1, SIZE(nhc % nvt, 1) 00627 ind = ind + 1 00628 nhc % nvt(j,i) % mass = buffer ( ind ) 00629 END DO 00630 END DO 00631 CALL section_vals_val_get(nose_section,"FORCE%_DEFAULT_KEYWORD_",r_vals=buffer,error=error) 00632 DO i = 1, SIZE ( nhc % nvt, 2) 00633 ind = map_info%index(i) 00634 ind = (ind-1) * nhc % nhc_len 00635 DO j = 1, SIZE(nhc % nvt, 1) 00636 ind = ind + 1 00637 nhc % nvt(j,i) % f = buffer ( ind ) 00638 END DO 00639 END DO 00640 END IF 00641 00642 IF (save_mem) THEN 00643 NULLIFY(work_section) 00644 work_section => section_vals_get_subs_vals(nose_section,"COORD",error=error) 00645 CALL section_vals_remove_values(work_section, error) 00646 NULLIFY(work_section) 00647 work_section => section_vals_get_subs_vals(nose_section,"VELOCITY",error=error) 00648 CALL section_vals_remove_values(work_section, error) 00649 NULLIFY(work_section) 00650 work_section => section_vals_get_subs_vals(nose_section,"FORCE",error=error) 00651 CALL section_vals_remove_values(work_section, error) 00652 NULLIFY(work_section) 00653 work_section => section_vals_get_subs_vals(nose_section,"MASS",error=error) 00654 CALL section_vals_remove_values(work_section, error) 00655 END IF 00656 00657 END IF 00658 00659 CALL timestop(handle) 00660 00661 END SUBROUTINE restart_nose 00662 00663 ! ***************************************************************************** 00669 SUBROUTINE init_nhc_variables ( nhc,temp_ext, para_env, globenv, error ) 00670 TYPE(lnhc_parameters_type), POINTER :: nhc 00671 REAL(KIND=dp), INTENT(IN) :: temp_ext 00672 TYPE(cp_para_env_type), POINTER :: para_env 00673 TYPE(global_environment_type), POINTER :: globenv 00674 TYPE(cp_error_type), INTENT(inout) :: error 00675 00676 CHARACTER(len=*), PARAMETER :: routineN = 'init_nhc_variables', 00677 routineP = moduleN//':'//routineN 00678 00679 INTEGER :: handle, i, icount, j, number, 00680 stat, tot_rn 00681 LOGICAL :: failure 00682 REAL(KIND=dp) :: akin, dum, temp 00683 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: array_of_rn 00684 TYPE(map_info_type), POINTER :: map_info 00685 00686 CALL timeset(routineN,handle) 00687 00688 failure = .FALSE. 00689 tot_rn = 0 00690 00691 ! first initializing the mass of the nhc variables 00692 nhc % nvt ( :, : ) % mass = nhc % nvt ( :, : ) % nkt * nhc%tau_nhc**2 00693 nhc % nvt ( :, : ) % eta = 0._dp 00694 nhc % nvt ( :, : ) % v = 0._dp 00695 nhc % nvt ( :, : ) % f = 0._dp 00696 00697 map_info => nhc%map_info 00698 SELECT CASE (map_info%dis_type) 00699 CASE ( do_thermo_only_master ) ! for NPT 00700 CASE DEFAULT 00701 tot_rn = nhc%glob_num_nhc*nhc%nhc_len 00702 00703 ALLOCATE(array_of_rn( tot_rn), STAT = stat) 00704 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00705 array_of_rn(:) = 0.0_dp 00706 END SELECT 00707 00708 SELECT CASE (map_info%dis_type ) 00709 CASE ( do_thermo_only_master ) ! for NPT 00710 ! Map deterministically determined random number to nhc % v 00711 DO i = 1, nhc%loc_num_nhc 00712 DO j = 1, nhc % nhc_len 00713 nhc%nvt(j,i)%v = next_random_number(globenv%gaussian_rng_stream,error=error) 00714 END DO 00715 END DO 00716 00717 akin = 0.0_dp 00718 DO i = 1, nhc%loc_num_nhc 00719 DO j = 1, nhc % nhc_len 00720 akin = akin + 0.5_dp * ( nhc % nvt ( j , i ) % mass * & 00721 nhc % nvt ( j , i ) % v * & 00722 nhc % nvt ( j , i ) % v ) 00723 END DO 00724 END DO 00725 number = nhc%loc_num_nhc 00726 00727 ! scale velocities to get the correct initial temperature 00728 temp = 2.0_dp*akin/REAL(number) 00729 IF (temp > 0.0_dp) temp = SQRT(temp_ext/temp) 00730 DO i = 1, nhc%loc_num_nhc 00731 DO j = 1, nhc % nhc_len 00732 nhc % nvt(j,i) % v = temp * nhc % nvt(j,i) % v 00733 nhc % nvt(j,i) % eta = 0.0_dp 00734 END DO 00735 END DO 00736 00737 ! initializing all of the forces on the thermostats 00738 DO i = 1, nhc%loc_num_nhc 00739 DO j = 2, nhc % nhc_len 00740 nhc % nvt(j,i) % f = nhc % nvt(j-1,i) % mass*nhc % nvt(j-1,i) % v* & 00741 nhc % nvt(j-1,i) % v - nhc % nvt(j,i) % nkt 00742 IF (nhc % nvt(j,i) % mass > 0.0_dp) THEN 00743 nhc % nvt(j,i) % f = nhc % nvt(j,i) % f/nhc % nvt(j,i) % mass 00744 END IF 00745 END DO 00746 END DO 00747 00748 CASE DEFAULT 00749 DO i=1,tot_rn 00750 array_of_rn(i) = next_random_number(globenv%gaussian_rng_stream,error=error) 00751 END DO 00752 ! Map deterministically determined random number to nhc % v 00753 DO i = 1, nhc%loc_num_nhc 00754 icount = map_info%index(i) 00755 icount = (icount-1) * nhc % nhc_len 00756 DO j = 1, nhc % nhc_len 00757 icount = icount + 1 00758 nhc % nvt(j,i) % v = array_of_rn ( icount ) 00759 ! WRITE ( *, * ) 'VEL', para_env%mepos, i,j, nhc%nvt(j,i)%v 00760 nhc % nvt(j,i) % eta = 0.0_dp 00761 END DO 00762 END DO 00763 DEALLOCATE ( array_of_rn, STAT = stat ) 00764 CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure) 00765 00766 number = nhc%glob_num_nhc 00767 CALL get_nhc_energies(nhc, dum, akin, para_env, error=error) 00768 00769 ! scale velocities to get the correct initial temperature 00770 temp = 2.0_dp*akin/REAL(number) 00771 IF (temp > 0.0_dp) temp = SQRT(temp_ext/temp) 00772 DO i = 1, nhc%loc_num_nhc 00773 DO j = 1, nhc % nhc_len 00774 nhc % nvt(j,i) % v = temp * nhc % nvt(j,i) % v 00775 END DO 00776 END DO 00777 00778 ! initializing all of the forces on the thermostats 00779 DO i = 1, nhc%loc_num_nhc 00780 DO j = 2, nhc % nhc_len 00781 nhc % nvt(j,i) % f = nhc % nvt(j-1,i) % mass*nhc % nvt(j-1,i) % v* & 00782 nhc % nvt(j-1,i) % v - nhc % nvt(j,i) % nkt 00783 IF (nhc % nvt(j,i) % mass > 0.0_dp) THEN 00784 nhc % nvt(j,i) % f = nhc % nvt(j,i) % f/nhc % nvt(j,i) % mass 00785 END IF 00786 END DO 00787 END DO 00788 00789 END SELECT 00790 00791 CALL timestop(handle) 00792 00793 END SUBROUTINE init_nhc_variables 00794 00795 ! ***************************************************************************** 00801 SUBROUTINE init_barostat_variables ( npt, tau_cell, temp_ext, nfree, ensemble, & 00802 cmass, globenv, error ) 00803 00804 TYPE(npt_info_type), DIMENSION(:, :), 00805 INTENT(INOUT) :: npt 00806 REAL(KIND=dp), INTENT(IN) :: tau_cell, temp_ext 00807 INTEGER, INTENT(IN) :: nfree, ensemble 00808 REAL(KIND=dp), INTENT(IN) :: cmass 00809 TYPE(global_environment_type), POINTER :: globenv 00810 TYPE(cp_error_type), INTENT(inout) :: error 00811 00812 CHARACTER(len=*), PARAMETER :: routineN = 'init_barostat_variables', 00813 routineP = moduleN//':'//routineN 00814 00815 INTEGER :: handle, i, j, number 00816 LOGICAL :: failure 00817 REAL(KIND=dp) :: akin, temp, v 00818 00819 CALL timeset(routineN,handle) 00820 00821 failure = .FALSE. 00822 temp = 0.0_dp 00823 00824 ! first initializing the mass of the nhc variables 00825 00826 npt(:,:)%eps = 0.0_dp 00827 npt(:,:)%v = 0.0_dp 00828 npt(:,:)%f = 0.0_dp 00829 SELECT CASE ( ensemble ) 00830 CASE ( npt_i_ensemble ) 00831 npt ( :, : ) % mass = REAL ( nfree + 3,KIND=dp) * temp_ext * tau_cell ** 2 00832 CASE ( npt_f_ensemble) 00833 npt ( :, : ) % mass = REAL ( nfree + 3,KIND=dp) * temp_ext * tau_cell ** 2 / 3.0_dp 00834 CASE ( nph_uniaxial_ensemble,nph_uniaxial_damped_ensemble ) 00835 npt ( :, : ) % mass = cmass 00836 CASE ( npe_f_ensemble) 00837 npt ( :, : ) % mass = REAL ( nfree + 3,KIND=dp) * temp_ext * tau_cell ** 2 / 3.0_dp 00838 CASE ( npe_i_ensemble ) 00839 npt ( :, : ) % mass = REAL ( nfree + 3,KIND=dp) * temp_ext * tau_cell ** 2 00840 END SELECT 00841 ! initializing velocities 00842 DO i = 1, SIZE ( npt,1) 00843 DO j = i, SIZE ( npt,2) 00844 v = next_random_number(globenv%gaussian_rng_stream,error=error) 00845 ! Symmetrizing the initial barostat velocities to ensure 00846 ! no rotation of the cell under NPT_F 00847 npt(j,i) % v = v 00848 npt(i,j) % v = v 00849 END DO 00850 END DO 00851 00852 akin = 0.0_dp 00853 DO i = 1, SIZE( npt,1 ) 00854 DO j = 1, SIZE( npt,2 ) 00855 akin = akin + 0.5_dp*( npt(j,i) % mass* npt(j,i)%v * npt(j,i)%v) 00856 END DO 00857 END DO 00858 00859 number = SIZE ( npt, 1 ) * SIZE ( npt, 2 ) 00860 00861 ! scale velocities to get the correct initial temperature 00862 IF ( number /= 0 ) THEN 00863 temp = 2.0_dp * akin / REAL ( number,KIND=dp) 00864 IF (temp > 0.0_dp) temp = SQRT(temp_ext/temp) 00865 ENDIF 00866 DO i = 1, SIZE(npt,1) 00867 DO j = i, SIZE(npt,2) 00868 npt(j,i) % v = temp * npt(j,i) % v 00869 npt(i,j) % v = npt(j,i) % v 00870 IF (debug_isotropic_limit) THEN 00871 npt(j,i) % v = 0.0_dp 00872 npt(i,j) % v = 0.0_dp 00873 WRITE (*, *) 'DEBUG ISOTROPIC LIMIT| INTIAL v_eps', npt(j,i) % v 00874 END IF 00875 END DO 00876 END DO 00877 00878 ! Zero barostat velocities for nph_uniaxial 00879 SELECT CASE ( ensemble ) 00880 ! Zero barostat velocities for nph_uniaxial 00881 CASE ( nph_uniaxial_ensemble, nph_uniaxial_damped_ensemble ) 00882 npt(:,:) % v = 0.0_dp 00883 END SELECT 00884 00885 CALL timestop(handle) 00886 00887 END SUBROUTINE init_barostat_variables 00888 00889 ! ***************************************************************************** 00893 SUBROUTINE init_nhc_forces ( nhc, error) 00894 00895 TYPE(lnhc_parameters_type), POINTER :: nhc 00896 TYPE(cp_error_type), INTENT(inout) :: error 00897 00898 CHARACTER(len=*), PARAMETER :: routineN = 'init_nhc_forces', 00899 routineP = moduleN//':'//routineN 00900 00901 INTEGER :: handle, i, j 00902 LOGICAL :: failure 00903 00904 CALL timeset(routineN,handle) 00905 00906 failure = .FALSE. 00907 00908 CPPrecondition(ASSOCIATED(nhc),cp_fatal_level,routineP,error,failure) 00909 ! assign the forces 00910 DO i = 1, SIZE ( nhc % nvt, 2 ) 00911 DO j = 2, SIZE ( nhc % nvt, 1 ) 00912 nhc % nvt ( j, i ) % f = nhc % nvt ( j-1, i ) % mass * & 00913 nhc % nvt ( j-1, i ) % v **2 - & 00914 nhc % nvt ( j, i ) % nkt 00915 IF (nhc % nvt ( j, i ) % mass > 0.0_dp) THEN 00916 nhc % nvt ( j, i ) % f = nhc % nvt ( j, i ) % f / nhc % nvt ( j, i ) % mass 00917 END IF 00918 END DO 00919 END DO 00920 00921 CALL timestop(handle) 00922 00923 END SUBROUTINE init_nhc_forces 00924 00925 END MODULE extended_system_init
1.7.3