|
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 hfx_types 00013 USE atomic_kind_types, ONLY: atomic_kind_type,& 00014 get_atomic_kind,& 00015 get_atomic_kind_set 00016 USE basis_set_types, ONLY: get_gto_basis_set,& 00017 gto_basis_set_type 00018 USE bibliography, ONLY: cite_reference,& 00019 guidon2008,& 00020 guidon2009 00021 USE cell_types, ONLY: cell_type,& 00022 get_cell,& 00023 plane_distance,& 00024 scaled_to_real 00025 USE cp_control_types, ONLY: dft_control_type 00026 USE cp_files, ONLY: close_file,& 00027 open_file 00028 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 00029 cp_print_key_unit_nr 00030 USE cp_para_types, ONLY: cp_para_env_type 00031 USE cp_units, ONLY: cp_unit_from_cp2k 00032 USE f77_blas 00033 USE hfx_helpers, ONLY: count_cells_perd,& 00034 next_image_cell_perd 00035 USE hfx_libint_wrapper, ONLY: initialize_libderiv,& 00036 initialize_libint,& 00037 terminate_libderiv,& 00038 terminate_libint 00039 USE hfx_libint_wrapper_types, ONLY: lib_deriv,& 00040 lib_int,& 00041 prim_data_f_size 00042 USE input_constants, ONLY: & 00043 do_hfx_auto_shells, do_hfx_potential_coulomb, & 00044 do_hfx_potential_gaussian, do_hfx_potential_id, do_hfx_potential_long, & 00045 do_hfx_potential_mix_cl, do_hfx_potential_mix_cl_trunc, & 00046 do_hfx_potential_mix_lg, do_hfx_potential_short, & 00047 do_hfx_potential_truncated, hfx_do_eval_energy, hfx_do_eval_forces, & 00048 use_aux_fit_basis_set 00049 USE input_section_types, ONLY: section_vals_get,& 00050 section_vals_get_subs_vals,& 00051 section_vals_type,& 00052 section_vals_val_get 00053 USE kinds, ONLY: default_path_length,& 00054 dp,& 00055 int_8 00056 USE machine, ONLY: m_chdir,& 00057 m_getcwd 00058 USE mathconstants, ONLY: pi 00059 USE orbital_pointers 00060 USE string_utilities, ONLY: compress 00061 USE t_c_g0, ONLY: free 00062 USE timings, ONLY: timeset,& 00063 timestop 00064 #include "cp_common_uses.h" 00065 00066 IMPLICIT NONE 00067 PRIVATE 00068 PUBLIC hfx_type, hfx_create, hfx_release,& 00069 hfx_set_distr_energy, & 00070 hfx_set_distr_forces, & 00071 hfx_cell_type, hfx_distribution, & 00072 hfx_potential_type, hfx_screening_type, hfx_periodic_type,& 00073 hfx_memory_type, hfx_load_balance_type, hfx_general_type,& 00074 hfx_compression_type,& 00075 hfx_container_type, hfx_container_node, hfx_cache_type, & 00076 hfx_basis_type, parse_memory_section, & 00077 hfx_init_container,& 00078 hfx_basis_info_type, hfx_screen_coeff_type,& 00079 hfx_reset_memory_usage_counter, pair_list_type, pair_list_element_type, & 00080 pair_set_list_type, hfx_p_kind, hfx_2D_map, hfx_pgf_list, & 00081 hfx_pgf_image, hfx_pgf_product_list, hfx_block_range_type, & 00082 alloc_containers, dealloc_containers, hfx_task_list_type, init_t_c_g0_lmax, & 00083 hfx_create_neighbor_cells, hfx_create_basis_types, hfx_release_basis_types 00084 00085 #define CACHE_SIZE 1024 00086 #define BITS_MAX_VAL 6 00087 00088 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hfx_types' 00089 INTEGER, PARAMETER, PUBLIC :: max_atom_block=32 00090 INTEGER, PARAMETER, PUBLIC :: max_images=27 00091 REAL(dp), PARAMETER, PUBLIC :: log_zero = -1000.0_dp 00092 REAL(dp), PARAMETER, PUBLIC :: powell_min_log = -20.0_dp 00093 REAL(KIND=dp), DIMENSION(0:10), 00094 PARAMETER, PUBLIC :: mul_fact = (/1.0_dp, 00095 1.1781_dp, 00096 1.3333_dp, 00097 1.4726_dp, 00098 1.6000_dp, 00099 1.7181_dp, 00100 1.8286_dp, 00101 1.9328_dp, 00102 2.0317_dp, 00103 2.1261_dp, 00104 2.2165_dp/) 00105 00106 INTEGER, SAVE :: init_t_c_g0_lmax = -1 00107 00108 !*** 00109 00110 ! ***************************************************************************** 00111 TYPE hfx_potential_type 00112 INTEGER :: potential_type !! 1/r/ erfc(wr)/r ... 00113 REAL(dp) :: omega !! w 00114 REAL(dp) :: scale_coulomb !! scaling factor for mixed potential 00115 REAL(dp) :: scale_longrange !! scaling factor for mixed potential 00116 REAL(dp) :: scale_gaussian !! scaling factor for mixed potential 00117 REAL(dp) :: cutoff_radius !! cutoff radius if cutoff potential in use 00118 CHARACTER(default_path_length) :: filename 00119 END TYPE 00120 00121 ! ***************************************************************************** 00122 TYPE hfx_screening_type 00123 REAL(dp) :: eps_schwarz !! threshold 00124 REAL(dp) :: eps_schwarz_forces !! threshold 00125 LOGICAL :: do_p_screening_forces !! screen on P^2 ? 00126 LOGICAL :: do_initial_p_screening !! screen on inital guess? 00127 END TYPE 00128 00129 ! ***************************************************************************** 00130 TYPE hfx_memory_type 00131 INTEGER :: max_memory !! user def max memory MB 00132 INTEGER(int_8) :: max_compression_counter !! corresponding number of reals 00133 INTEGER(int_8) :: final_comp_counter_energy 00134 LOGICAL :: do_all_on_the_fly !! max mem == 0 ? 00135 REAL(dp) :: eps_storage_scaling 00136 INTEGER :: cache_size 00137 INTEGER :: bits_max_val 00138 INTEGER :: actual_memory_usage 00139 INTEGER :: actual_memory_usage_disk 00140 INTEGER(int_8) :: max_compression_counter_disk 00141 LOGICAL :: do_disk_storage 00142 CHARACTER(len=default_path_length) :: storage_location 00143 INTEGER(int_8) :: ram_counter 00144 INTEGER(int_8) :: ram_counter_forces 00145 INTEGER(int_8) :: size_p_screen 00146 LOGICAL :: treat_forces_in_core 00147 LOGICAL :: recalc_forces 00148 END TYPE 00149 00150 ! ***************************************************************************** 00151 TYPE hfx_periodic_type 00152 INTEGER :: number_of_shells !! number of periodic image cells 00153 LOGICAL :: do_periodic !! periodic ? 00154 INTEGER :: perd(3) !! x,xy,xyz,... 00155 INTEGER :: mode 00156 REAL(dp) :: R_max_stress 00157 INTEGER :: number_of_shells_from_input 00158 END TYPE 00159 00160 ! ***************************************************************************** 00161 TYPE hfx_load_balance_type 00162 INTEGER :: nbins 00163 INTEGER :: block_size 00164 INTEGER :: nblocks 00165 LOGICAL :: rtp_redistribute 00166 LOGICAL :: blocks_initialized 00167 LOGICAL :: do_randomize 00168 END TYPE 00169 00170 ! ***************************************************************************** 00171 TYPE hfx_general_type 00172 REAL(dp) :: fraction !! for hybrids 00173 LOGICAL :: treat_lsd_in_core 00174 END TYPE 00175 00176 ! ***************************************************************************** 00177 TYPE hfx_cell_type 00178 REAL(dp) :: cell(3) 00179 REAL(dp) :: cell_r(3) 00180 END TYPE 00181 00182 ! ***************************************************************************** 00183 TYPE hfx_distribution 00184 INTEGER(int_8) :: istart 00185 INTEGER(int_8) :: number_of_atom_quartets 00186 INTEGER(int_8) :: cost 00187 REAL(KIND=dp) :: time_first_scf 00188 REAL(KIND=dp) :: time_other_scf 00189 REAL(KIND=dp) :: time_forces 00190 INTEGER(int_8) :: ram_counter 00191 END TYPE 00192 00193 ! ***************************************************************************** 00194 TYPE pair_list_element_type 00195 INTEGER, DIMENSION(2) :: pair 00196 INTEGER, DIMENSION(2) :: set_bounds 00197 INTEGER, DIMENSION(2) :: kind_pair 00198 REAL(KIND=dp) :: r1(3),r2(3) 00199 REAL(KIND=dp) :: dist2 00200 END TYPE 00201 00202 ! ***************************************************************************** 00203 TYPE pair_set_list_type 00204 INTEGER, DIMENSION(2) :: pair 00205 END TYPE 00206 00207 ! ***************************************************************************** 00208 TYPE pair_list_type 00209 TYPE(pair_list_element_type), DIMENSION(max_atom_block**2) :: elements 00210 INTEGER :: n_element 00211 END TYPE pair_list_type 00212 00213 ! ***************************************************************************** 00214 TYPE hfx_compression_type 00215 INTEGER(int_8), DIMENSION(:), POINTER :: compressed_data 00216 INTEGER(int_8) :: element_pointer 00217 INTEGER :: bit_pointer 00218 END TYPE 00219 00220 ! ***************************************************************************** 00221 TYPE hfx_cache_type 00222 INTEGER(int_8), DIMENSION(CACHE_SIZE) :: DATA 00223 INTEGER :: element_counter 00224 END TYPE 00225 00226 ! ***************************************************************************** 00227 TYPE hfx_container_node 00228 TYPE(hfx_container_node), POINTER :: next, prev 00229 INTEGER(int_8), DIMENSION(CACHE_SIZE) :: DATA 00230 END TYPE 00231 00232 ! ***************************************************************************** 00233 TYPE hfx_container_type 00234 TYPE(hfx_container_node), POINTER :: first, current 00235 INTEGER :: element_counter 00236 INTEGER(int_8) :: file_counter 00237 CHARACTER(LEN=5) :: desc 00238 INTEGER :: unit 00239 CHARACTER(default_path_length) :: filename 00240 END TYPE 00241 00242 ! ***************************************************************************** 00243 TYPE hfx_basis_type 00244 INTEGER, DIMENSION(:), POINTER :: lmax 00245 INTEGER, DIMENSION(:), POINTER :: lmin 00246 INTEGER, DIMENSION(:), POINTER :: npgf 00247 INTEGER :: nset 00248 REAL(dp), DIMENSION(:,:), POINTER :: zet 00249 INTEGER, DIMENSION(:), POINTER :: nsgf 00250 INTEGER, DIMENSION(:,:), POINTER :: first_sgf 00251 REAL(dp), DIMENSION(:,:), POINTER :: sphi 00252 INTEGER :: nsgf_total 00253 INTEGER, DIMENSION(:,:), POINTER :: nl 00254 INTEGER, DIMENSION(:,:), POINTER :: nsgfl 00255 INTEGER, DIMENSION(:), POINTER :: nshell 00256 REAL(dp), DIMENSION(:,:,:,:), POINTER 00257 :: sphi_ext 00258 REAL(dp), DIMENSION(:), POINTER :: set_radius 00259 REAL(dp), DIMENSION(:,:), POINTER :: pgf_radius 00260 REAL(dp) :: kind_radius 00261 END TYPE 00262 00263 ! ***************************************************************************** 00264 TYPE hfx_basis_info_type 00265 INTEGER :: max_set 00266 INTEGER :: max_sgf 00267 INTEGER :: max_am 00268 END TYPE 00269 00270 ! ***************************************************************************** 00271 TYPE hfx_screen_coeff_type 00272 REAL(dp) :: x(2) 00273 END TYPE 00274 00275 ! ***************************************************************************** 00276 TYPE hfx_p_kind 00277 REAL(dp), DIMENSION(:,:,:,:), POINTER :: p_kind 00278 END TYPE 00279 00280 ! ***************************************************************************** 00281 TYPE hfx_2D_map 00282 INTEGER, DIMENSION(:), POINTER :: iatom_list 00283 INTEGER, DIMENSION(:), POINTER :: jatom_list 00284 END TYPE 00285 00286 ! ***************************************************************************** 00287 TYPE hfx_pgf_image 00288 REAL(dp) :: ra(3), rb(3) 00289 REAL(dp) :: rab2 00290 REAL(dp) :: S1234 00291 REAL(dp) :: P(3) 00292 REAL(dp) :: R 00293 REAL(dp) :: pgf_max 00294 END TYPE 00295 00296 ! ***************************************************************************** 00297 TYPE hfx_pgf_list 00298 TYPE(hfx_pgf_image), DIMENSION(:), POINTER 00299 :: image_list 00300 INTEGER :: nimages 00301 REAL(dp) :: zetapzetb 00302 REAL(dp) :: ZetaInv 00303 REAL(dp) :: zeta, zetb 00304 INTEGER :: ipgf, jpgf 00305 END TYPE 00306 00307 ! ***************************************************************************** 00308 TYPE hfx_pgf_product_list 00309 REAL(dp) :: ra(3), rb(3), rc(3), rd(3) 00310 REAL(dp) :: ZetapEtaInv 00311 REAL(dp) :: Rho, RhoInv 00312 REAL(dp) :: P(3), Q(3), W(3) 00313 REAL(dp) :: AB(3), CD(3) 00314 REAL(dp) :: Fm(prim_data_f_size) 00315 END TYPE 00316 00317 ! ***************************************************************************** 00318 TYPE hfx_block_range_type 00319 INTEGER :: istart,iend 00320 INTEGER(int_8) :: cost 00321 END TYPE 00322 00323 ! ***************************************************************************** 00324 TYPE hfx_task_list_type 00325 INTEGER :: thread_id 00326 INTEGER :: bin_id 00327 INTEGER(int_8) :: cost 00328 END TYPE 00329 00330 ! ***************************************************************************** 00365 TYPE hfx_type 00366 TYPE(hfx_potential_type) :: potential_parameter 00367 TYPE(hfx_screening_type) :: screening_parameter 00368 TYPE(hfx_memory_type) :: memory_parameter 00369 TYPE(hfx_periodic_type) :: periodic_parameter 00370 TYPE(hfx_load_balance_type) :: load_balance_parameter 00371 TYPE(hfx_general_type) :: general_parameter 00372 TYPE(hfx_container_type), DIMENSION(:), 00373 POINTER :: maxval_container 00374 TYPE(hfx_cache_type), DIMENSION(:), 00375 POINTER :: maxval_cache 00376 TYPE(hfx_container_type), DIMENSION(:,:), 00377 POINTER :: integral_containers 00378 TYPE(hfx_cache_type), DIMENSION(:,:), 00379 POINTER :: integral_caches 00380 TYPE(hfx_container_type), DIMENSION(:), 00381 POINTER :: maxval_container_forces 00382 TYPE(hfx_cache_type), DIMENSION(:), 00383 POINTER :: maxval_cache_forces 00384 TYPE(hfx_container_type), DIMENSION(:,:), 00385 POINTER :: integral_containers_forces 00386 TYPE(hfx_cache_type), DIMENSION(:,:), 00387 POINTER :: integral_caches_forces 00388 00389 TYPE(hfx_container_type), POINTER :: maxval_container_disk 00390 TYPE(hfx_cache_type) :: maxval_cache_disk 00391 TYPE(hfx_container_type), POINTER, 00392 DIMENSION(:) :: integral_containers_disk 00393 TYPE(hfx_cache_type) :: integral_caches_disk(64) 00394 00395 TYPE(hfx_cell_type), DIMENSION(:), 00396 POINTER :: neighbor_cells 00397 TYPE(hfx_distribution), DIMENSION(:), 00398 POINTER :: distribution_energy 00399 TYPE(hfx_distribution), DIMENSION(:), 00400 POINTER :: distribution_forces 00401 INTEGER, DIMENSION(:,:), POINTER :: is_assoc_atomic_block 00402 INTEGER :: number_of_p_entries 00403 TYPE(hfx_basis_type), DIMENSION(:), 00404 POINTER :: basis_parameter 00405 INTEGER :: n_rep_hf 00406 LOGICAL :: b_first_load_balance_energy, 00407 b_first_load_balance_forces 00408 REAL(dp), DIMENSION(:), POINTER :: full_ks_alpha 00409 REAL(dp), DIMENSION(:), POINTER :: full_ks_beta 00410 TYPE(lib_int) :: lib 00411 TYPE(hfx_basis_info_type) :: basis_info 00412 TYPE(hfx_screen_coeff_type), 00413 DIMENSION(:,:,:,:,:,:), POINTER :: screen_funct_coeffs_pgf, 00414 pair_dist_radii_pgf 00415 TYPE(hfx_screen_coeff_type), 00416 DIMENSION(:,:,:,:), POINTER :: screen_funct_coeffs_set 00417 TYPE(hfx_screen_coeff_type), 00418 DIMENSION(:,:), POINTER :: screen_funct_coeffs_kind 00419 LOGICAL :: screen_funct_is_initialized 00420 TYPE(hfx_p_kind), DIMENSION(:), POINTER :: initial_p 00421 TYPE(hfx_p_kind), DIMENSION(:), POINTER :: initial_p_forces 00422 INTEGER, DIMENSION(:), POINTER :: map_atom_to_kind_atom 00423 TYPE(hfx_2D_map), DIMENSION(:), POINTER :: map_atoms_to_cpus 00424 INTEGER, DIMENSION(:,:), POINTER :: atomic_block_offset 00425 INTEGER, DIMENSION(:,:,:,:), POINTER :: set_offset 00426 INTEGER, DIMENSION(:), POINTER :: block_offset 00427 TYPE(hfx_block_range_type), DIMENSION(:), 00428 POINTER :: blocks 00429 TYPE(hfx_task_list_type), DIMENSION(:), 00430 POINTER :: task_list 00431 REAL(dp), DIMENSION(:,:), POINTER :: pmax_atom, pmax_atom_forces 00432 TYPE(lib_int) :: lib_int 00433 TYPE(lib_deriv) :: lib_deriv 00434 REAL(dp), DIMENSION(:,:), POINTER :: pmax_block 00435 LOGICAL, DIMENSION(:,:), POINTER :: atomic_pair_list 00436 END TYPE hfx_type 00437 00438 CONTAINS 00439 00440 ! ***************************************************************************** 00454 SUBROUTINE hfx_create(x_data,para_env, hfx_section,natom,atomic_kind_set, & 00455 dft_control, cell, error) 00456 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data 00457 TYPE(cp_para_env_type) :: para_env 00458 TYPE(section_vals_type), POINTER :: hfx_section 00459 INTEGER, INTENT(IN) :: natom 00460 TYPE(atomic_kind_type), DIMENSION(:), 00461 POINTER :: atomic_kind_set 00462 TYPE(dft_control_type), POINTER :: dft_control 00463 TYPE(cell_type), POINTER :: cell 00464 TYPE(cp_error_type), INTENT(inout) :: error 00465 00466 CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_create', 00467 routineP = moduleN//':'//routineN 00468 00469 CHARACTER(LEN=512) :: error_msg 00470 CHARACTER(LEN=default_path_length) :: char_val 00471 INTEGER :: handle, i, i_thread, iatom, ikind, int_val, irep, jkind, 00472 max_set, n_rep_hf, n_threads, natom_a, natom_b, nkind, nseta, nsetb, 00473 pbc_shells, stat, storage_id 00474 INTEGER, DIMENSION(:), POINTER :: atom2kind, kind_of 00475 LOGICAL :: exist, failure, logic_val 00476 REAL(dp) :: real_val 00477 TYPE(atomic_kind_type), POINTER :: atom_kind 00478 TYPE(hfx_type), POINTER :: actual_x_data 00479 TYPE(section_vals_type), POINTER :: hf_pbc_section, hf_sub_section 00480 00481 !$ INTEGER :: omp_get_max_threads 00482 00483 CALL timeset(routineN,handle) 00484 00485 failure = .FALSE. 00486 CALL cite_reference(Guidon2008) 00487 CALL cite_reference(Guidon2009) 00488 !! There might be 2 hf sections 00489 CALL section_vals_get(hfx_section,n_repetition=n_rep_hf,error=error) 00490 n_threads = 1 00491 !$ n_threads = omp_get_max_threads() 00492 ALLOCATE(x_data(n_rep_hf,n_threads),STAT=stat) 00493 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00494 DO i_thread=1,n_threads 00495 DO irep=1,n_rep_hf 00496 actual_x_data => x_data(irep,i_thread) 00497 !! Get data from input file 00498 !! 00499 !! GENERAL params 00500 CALL section_vals_val_get(hfx_section, "FRACTION", r_val=real_val, i_rep_section=irep,error=error) 00501 actual_x_data%general_parameter%fraction = real_val 00502 actual_x_data%n_rep_hf = n_rep_hf 00503 00504 NULLIFY(actual_x_data%map_atoms_to_cpus) 00505 00506 CALL section_vals_val_get(hfx_section, "TREAT_LSD_IN_CORE", l_val=logic_val, i_rep_section=irep,error=error) 00507 actual_x_data%general_parameter%treat_lsd_in_core = logic_val 00508 00509 !! MEMORY section 00510 hf_sub_section => section_vals_get_subs_vals(hfx_section,"MEMORY",i_rep_section=irep, error=error) 00511 CALL parse_memory_section(actual_x_data%memory_parameter, hf_sub_section, storage_id, i_thread,& 00512 n_threads, para_env, irep, skip_disk=.FALSE., skip_in_core_forces=.FALSE., error=error) 00513 00514 !! PERIODIC section 00515 hf_sub_section => section_vals_get_subs_vals(hfx_section,"PERIODIC", i_rep_section=irep, error=error) 00516 CALL section_vals_val_get(hf_sub_section, "NUMBER_OF_SHELLS", i_val=int_val, error=error) 00517 actual_x_data%periodic_parameter%number_of_shells = int_val 00518 actual_x_data%periodic_parameter%mode = int_val 00519 CALL get_cell(cell=cell,periodic=actual_x_data%periodic_parameter%perd) 00520 IF( SUM(actual_x_data%periodic_parameter%perd) == 0 ) THEN 00521 actual_x_data%periodic_parameter%do_periodic = .FALSE. 00522 ELSE 00523 actual_x_data%periodic_parameter%do_periodic = .TRUE. 00524 END IF 00525 00526 !! SCREENING section 00527 hf_sub_section => section_vals_get_subs_vals(hfx_section,"SCREENING",i_rep_section=irep, error=error) 00528 CALL section_vals_val_get(hf_sub_section, "EPS_SCHWARZ", r_val=real_val, error=error) 00529 actual_x_data%screening_parameter%eps_schwarz = real_val 00530 CALL section_vals_val_get(hf_sub_section, "EPS_SCHWARZ_FORCES", r_val=real_val, error=error) 00531 actual_x_data%screening_parameter%eps_schwarz_forces = real_val 00532 CALL section_vals_val_get(hf_sub_section, "SCREEN_P_FORCES", l_val=logic_val, error=error) 00533 actual_x_data%screening_parameter%do_p_screening_forces = logic_val 00534 CALL section_vals_val_get(hf_sub_section, "SCREEN_ON_INITIAL_P", l_val=logic_val, error=error) 00535 actual_x_data%screening_parameter%do_initial_p_screening = logic_val 00536 actual_x_data%screen_funct_is_initialized = .FALSE. 00537 00538 !! INTERACTION_POTENTIAL section 00539 hf_sub_section => section_vals_get_subs_vals(hfx_section,"INTERACTION_POTENTIAL",i_rep_section=irep,error=error) 00540 CALL section_vals_val_get(hf_sub_section, "POTENTIAL_TYPE", i_val=int_val, error=error) 00541 actual_x_data%potential_parameter%potential_type = int_val 00542 CALL section_vals_val_get(hf_sub_section, "OMEGA", r_val=real_val, error=error) 00543 actual_x_data%potential_parameter%omega = real_val 00544 CALL section_vals_val_get(hf_sub_section, "SCALE_COULOMB", r_val=real_val, error=error) 00545 actual_x_data%potential_parameter%scale_coulomb = real_val 00546 CALL section_vals_val_get(hf_sub_section, "SCALE_LONGRANGE", r_val=real_val, error=error) 00547 actual_x_data%potential_parameter%scale_longrange = real_val 00548 CALL section_vals_val_get(hf_sub_section, "SCALE_GAUSSIAN", r_val=real_val, error=error) 00549 actual_x_data%potential_parameter%scale_gaussian = real_val 00550 IF (actual_x_data%potential_parameter%potential_type == do_hfx_potential_truncated .OR. & 00551 actual_x_data%potential_parameter%potential_type == do_hfx_potential_mix_cl_trunc) THEN 00552 CALL section_vals_val_get(hf_sub_section, "CUTOFF_RADIUS", r_val=real_val, error=error) 00553 actual_x_data%potential_parameter%cutoff_radius = real_val 00554 CALL section_vals_val_get(hf_sub_section, "T_C_G_DATA", c_val=char_val, error=error) 00555 CALL compress(char_val,.TRUE.) 00556 ! ** Check if file is there 00557 INQUIRE(FILE=char_val,exist=exist) 00558 IF ( .NOT. exist ) THEN 00559 WRITE(error_msg,'(A,A,A)') "Truncated hfx calculation requested. The file containing "//& 00560 "the data could not be found at ", TRIM(char_val), " Please check T_C_G_DATA "//& 00561 "in the INTERCATION_POTENTIAL section" 00562 CALL cp_assert( .FALSE. , cp_failure_level,cp_assertion_failed,routineP,& 00563 error_msg , & 00564 error,& 00565 only_ionode=.TRUE.) 00566 ELSE 00567 actual_x_data%potential_parameter%filename = char_val 00568 END IF 00569 END IF 00570 IF (actual_x_data%potential_parameter%potential_type == do_hfx_potential_short) THEN 00571 CALL erfc_cutoff(actual_x_data%screening_parameter%eps_schwarz,& 00572 actual_x_data%potential_parameter%omega,& 00573 actual_x_data%potential_parameter%cutoff_radius, error) 00574 END IF 00575 00576 !! LOAD_BALANCE section 00577 hf_sub_section => section_vals_get_subs_vals(hfx_section,"LOAD_BALANCE",i_rep_section=irep, error=error) 00578 CALL section_vals_val_get(hf_sub_section, "NBINS", i_val=int_val, error=error) 00579 actual_x_data%load_balance_parameter%nbins = MAX(1,int_val) 00580 actual_x_data%load_balance_parameter%blocks_initialized = .FALSE. 00581 00582 CALL section_vals_val_get(hf_sub_section, "RANDOMIZE", l_val=logic_val, error=error) 00583 actual_x_data%load_balance_parameter%do_randomize = logic_val 00584 00585 actual_x_data%load_balance_parameter%rtp_redistribute=.FALSE. 00586 IF(ASSOCIATED(dft_control%rtp_control))& 00587 actual_x_data%load_balance_parameter%rtp_redistribute=dft_control%rtp_control%hfx_redistribute 00588 00589 CALL section_vals_val_get(hf_sub_section, "BLOCK_SIZE", i_val=int_val, error=error) 00590 ! negative values ask for a computed default 00591 IF (int_val<=0) THEN 00592 ! this gives a reasonable number of blocks for binning, yet typically results in blocking. 00593 int_val = CEILING( 0.1_dp * natom / & 00594 REAL(actual_x_data%load_balance_parameter%nbins * n_threads * para_env%num_pe ,KIND=dp)**(0.25_dp)) 00595 ENDIF 00596 ! at least 1 atom per block, and avoid overly large blocks 00597 actual_x_data%load_balance_parameter%block_size = MIN(max_atom_block, MAX(1, int_val)) 00598 00599 CALL hfx_create_basis_types(actual_x_data%basis_parameter,actual_x_data%basis_info, atomic_kind_set,& 00600 dft_control%do_admm, error) 00601 00602 !!************************************************************************************************** 00603 !! ** !! ** This code writes the contraction routines 00604 !! ** !! ** Very UGLY: BASIS_SET has to be 1 primitive and lmin=lmax=l. For g-functions 00605 !! ** !! ** 00606 !! ** !! ** 1 4 4 1 1 00607 !! ** !! ** 1.0 1.0 00608 !! ** !! ** 00609 !! ** k = max_am - 1 00610 !! ** write(filename,'(A,I0,A)') "sphi",k+1,"a" 00611 !! ** OPEN(UNIT=31415,FILE=filename) 00612 !! ** DO i=ncoset(k)+1,SIZE(sphi_a,1) 00613 !! ** DO j=1,SIZE(sphi_a,2) 00614 !! ** IF( sphi_a(i,j) /= 0.0_dp) THEN 00615 !! ** write(31415,'(A,I0,A,I0,A,I0,A,I0,A,I0,A)') "buffer1(i+imax*(",& 00616 !! ** j,& 00617 !! ** "-1)) = buffer1(i+imax*(",& 00618 !! ** j,& 00619 !! ** "-1)) + work(",& 00620 !! ** i-ncoset(k),& 00621 !! ** "+(i-1)*kmax) * sphi_a(",& 00622 !! ** i-ncoset(k),& 00623 !! ** ",",& 00624 !! ** j,& 00625 !! ** "+s_offset_a1)" 00626 !! ** END IF 00627 !! ** END DO 00628 !! ** END DO 00629 !! ** CLOSE(UNIT=31415) 00630 !! ** write(filename,'(A,I0,A)') "sphi",k+1,"b" 00631 !! ** OPEN(UNIT=31415,FILE=filename) 00632 !! ** DO i=ncoset(k)+1,SIZE(sphi_a,1) 00633 !! ** DO j=1,SIZE(sphi_a,2) 00634 !! ** IF( sphi_a(i,j) /= 0.0_dp) THEN 00635 !! ** write(31415,'(A,I0,A,I0,A,I0,A,I0,A,I0,A)') "buffer2(i+imax*(",& 00636 !! ** j,& 00637 !! ** "-1)) = buffer2(i+imax*(",& 00638 !! ** j,& 00639 !! ** "-1)) + buffer1(",& 00640 !! ** i-ncoset(k),& 00641 !! ** "+(i-1)*kmax) * sphi_b(",& 00642 !! ** i-ncoset(k),& 00643 !! ** ",",& 00644 !! ** j,& 00645 !! ** "+s_offset_b1)" 00646 !! ** 00647 !! ** END IF 00648 !! ** END DO 00649 !! ** END DO 00650 !! ** CLOSE(UNIT=31415) 00651 !! ** write(filename,'(A,I0,A)') "sphi",k+1,"c" 00652 !! ** OPEN(UNIT=31415,FILE=filename) 00653 !! ** DO i=ncoset(k)+1,SIZE(sphi_a,1) 00654 !! ** DO j=1,SIZE(sphi_a,2) 00655 !! ** IF( sphi_a(i,j) /= 0.0_dp) THEN 00656 !! ** write(31415,'(A,I0,A,I0,A,I0,A,I0,A,I0,A)') "buffer1(i+imax*(",& 00657 !! ** j,& 00658 !! ** "-1)) = buffer1(i+imax*(",& 00659 !! ** j,& 00660 !! ** "-1)) + buffer2(",& 00661 !! ** i-ncoset(k),& 00662 !! ** "+(i-1)*kmax) * sphi_c(",& 00663 !! ** i-ncoset(k),& 00664 !! ** ",",& 00665 !! ** j,& 00666 !! ** "+s_offset_c1)" 00667 !! ** 00668 !! ** END IF 00669 !! ** END DO 00670 !! ** END DO 00671 !! ** CLOSE(UNIT=31415) 00672 !! ** write(filename,'(A,I0,A)') "sphi",k+1,"d" 00673 !! ** OPEN(UNIT=31415,FILE=filename) 00674 !! ** DO i=ncoset(k)+1,SIZE(sphi_a,1) 00675 !! ** DO j=1,SIZE(sphi_a,2) 00676 !! ** IF( sphi_a(i,j) /= 0.0_dp) THEN 00677 !! ** 00678 !! ** 00679 !! ** write(31415,'(A,I0,A)') "primitives(s_offset_a1+i3, s_offset_b1+i2, s_offset_c1+i1, s_offset_d1+",& 00680 !! ** j,")= &" 00681 !! ** write(31415,'(A,I0,A)') "primitives(s_offset_a1+i3, s_offset_b1+i2, s_offset_c1+i1, s_offset_d1+",& 00682 !! ** j,")+ &" 00683 !! ** write(31415,'(A,I0,A,I0,A,I0,A)') "buffer1(",& 00684 !! ** i-ncoset(k),& 00685 !! ** "+(i-1)*kmax) * sphi_d(",& 00686 !! ** i-ncoset(k),& 00687 !! ** ",",& 00688 !! ** j,& 00689 !! ** "+s_offset_d1)" 00690 !! ** 00691 !! ** 00692 !! ** END IF 00693 !! ** END DO 00694 !! ** END DO 00695 !! ** CLOSE(UNIT=31415) 00696 !! ** stop 00697 !! **************************************************************************************************** 00698 00699 IF(actual_x_data%periodic_parameter%do_periodic) THEN 00700 hf_pbc_section => section_vals_get_subs_vals(hfx_section,"PERIODIC",i_rep_section=irep, error=error) 00701 CALL section_vals_val_get(hf_pbc_section,"NUMBER_OF_SHELLS",i_val=pbc_shells,error=error) 00702 actual_x_data%periodic_parameter%number_of_shells_from_input = pbc_shells 00703 ALLOCATE(actual_x_data%neighbor_cells(1),STAT=stat) 00704 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00705 CALL hfx_create_neighbor_cells(actual_x_data,pbc_shells, cell, i_thread, error=error) 00706 ELSE 00707 ALLOCATE(actual_x_data%neighbor_cells(1),STAT=stat) 00708 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00709 ! ** Initialize this guy to enable non periodic stress regtests 00710 actual_x_data%periodic_parameter%R_max_stress = 1.0_dp 00711 END IF 00712 00713 nkind = SIZE(atomic_kind_set,1) 00714 max_set=actual_x_data%basis_info%max_set 00715 00716 !! ** This guy is allocated on the master thread only 00717 IF (i_thread == 1 ) THEN 00718 ALLOCATE(actual_x_data%is_assoc_atomic_block(natom,natom),STAT=stat) 00719 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00720 ALLOCATE(actual_x_data%atomic_block_offset(natom,natom),STAT=stat) 00721 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00722 ALLOCATE(actual_x_data%set_offset(max_set, max_set, nkind, nkind),STAT=stat) 00723 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00724 ALLOCATE(actual_x_data%block_offset(para_env%num_pe+1),STAT=stat) 00725 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00726 END IF 00727 00728 ALLOCATE(actual_x_data%distribution_forces(1),STAT=stat) 00729 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00730 ALLOCATE(actual_x_data%distribution_energy(1),STAT=stat) 00731 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00732 00733 actual_x_data%memory_parameter%size_p_screen = 0_int_8 00734 IF( i_thread == 1 ) THEN 00735 ALLOCATE(actual_x_data%atomic_pair_list(natom,natom),STAT=stat) 00736 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00737 END IF 00738 00739 IF(actual_x_data%screening_parameter%do_initial_p_screening .OR. & 00740 actual_x_data%screening_parameter%do_p_screening_forces) THEN 00741 !! ** This guy is allocated on the master thread only 00742 IF (i_thread == 1 ) THEN 00743 ALLOCATE(actual_x_data%pmax_atom(natom,natom),STAT=stat) 00744 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00745 ALLOCATE(actual_x_data%initial_p(nkind*(nkind+1)/2),STAT=stat) 00746 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00747 i = 1 00748 DO ikind=1,nkind 00749 NULLIFY(atom_kind) 00750 atom_kind => atomic_kind_set(ikind) 00751 CALL get_atomic_kind(atomic_kind=atom_kind,& 00752 natom=natom_a) 00753 nseta=actual_x_data%basis_parameter(ikind)%nset 00754 DO jkind=ikind,nkind 00755 NULLIFY(atom_kind) 00756 atom_kind => atomic_kind_set(jkind) 00757 CALL get_atomic_kind(atomic_kind=atom_kind,& 00758 natom=natom_b) 00759 nsetb=actual_x_data%basis_parameter(jkind)%nset 00760 ALLOCATE(actual_x_data%initial_p(i)%p_kind(nseta, nsetb, natom_a, natom_b),STAT=stat) 00761 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00762 actual_x_data%memory_parameter%size_p_screen = & 00763 actual_x_data%memory_parameter%size_p_screen + nseta*nsetb*natom_a*natom_b 00764 i = i + 1 00765 END DO 00766 END DO 00767 IF( actual_x_data%memory_parameter%treat_forces_in_core ) THEN 00768 ALLOCATE(actual_x_data%pmax_atom_forces(natom,natom),STAT=stat) 00769 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00770 00771 ALLOCATE(actual_x_data%initial_p_forces(nkind*(nkind+1)/2),STAT=stat) 00772 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00773 i = 1 00774 DO ikind=1,nkind 00775 NULLIFY(atom_kind) 00776 atom_kind => atomic_kind_set(ikind) 00777 CALL get_atomic_kind(atomic_kind=atom_kind,& 00778 natom=natom_a) 00779 nseta=actual_x_data%basis_parameter(ikind)%nset 00780 DO jkind=ikind,nkind 00781 NULLIFY(atom_kind) 00782 atom_kind => atomic_kind_set(jkind) 00783 CALL get_atomic_kind(atomic_kind=atom_kind,& 00784 natom=natom_b) 00785 nsetb=actual_x_data%basis_parameter(jkind)%nset 00786 ALLOCATE(actual_x_data%initial_p_forces(i)%p_kind(nseta, nsetb, natom_a, natom_b),STAT=stat) 00787 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00788 actual_x_data%memory_parameter%size_p_screen = & 00789 actual_x_data%memory_parameter%size_p_screen + nseta*nsetb*natom_a*natom_b 00790 i = i + 1 00791 END DO 00792 END DO 00793 END IF 00794 END IF 00795 ALLOCATE(actual_x_data%map_atom_to_kind_atom(natom),STAT=stat) 00796 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00797 ALLOCATE(kind_of(natom),STAT=stat) 00798 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00799 ALLOCATE (atom2kind(nkind),STAT=stat) 00800 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00801 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 00802 kind_of=kind_of) 00803 00804 atom2kind = 0 00805 DO iatom=1,natom 00806 ikind = kind_of(iatom) 00807 atom2kind(ikind) = atom2kind(ikind) + 1 00808 actual_x_data%map_atom_to_kind_atom(iatom) = atom2kind(ikind) 00809 END DO 00810 DEALLOCATE(kind_of, atom2kind, STAT=stat) 00811 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00812 END IF 00813 00814 ! ** Initialize libint type 00815 CALL initialize_libint(actual_x_data%lib, actual_x_data%basis_info%max_am, error) 00816 CALL initialize_libderiv(actual_x_data%lib_deriv, actual_x_data%basis_info%max_am, error) 00817 00818 CALL alloc_containers(actual_x_data, 1, hfx_do_eval_energy, error) 00819 CALL alloc_containers(actual_x_data, 1, hfx_do_eval_forces, error) 00820 00821 actual_x_data%maxval_cache_disk%element_counter = 1 00822 ALLOCATE(actual_x_data%maxval_container_disk,STAT=stat) 00823 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00824 ALLOCATE(actual_x_data%maxval_container_disk%first,STAT=stat) 00825 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00826 actual_x_data%maxval_container_disk%first%prev => NULL() 00827 actual_x_data%maxval_container_disk%first%next => NULL() 00828 actual_x_data%maxval_container_disk%current => actual_x_data%maxval_container_disk%first 00829 actual_x_data%maxval_container_disk%current%data = 0 00830 actual_x_data%maxval_container_disk%element_counter = 1 00831 actual_x_data%maxval_container_disk%file_counter = 1 00832 actual_x_data%maxval_container_disk%desc = 'Max_' 00833 actual_x_data%maxval_container_disk%unit = -1 00834 WRITE(actual_x_data%maxval_container_disk%filename,'(A,I0,A,A,A)') TRIM(actual_x_data%memory_parameter%storage_location), & 00835 storage_id,"_",actual_x_data%maxval_container_disk%desc, "6" 00836 CALL compress(actual_x_data%maxval_container_disk%filename, .TRUE.) 00837 ALLOCATE(actual_x_data%integral_containers_disk(64), STAT=stat) 00838 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00839 DO i=1,64 00840 actual_x_data%integral_caches_disk(i)%element_counter = 1 00841 actual_x_data%integral_caches_disk(i)%data = 0 00842 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00843 ALLOCATE(actual_x_data%integral_containers_disk(i)%first,STAT=stat) 00844 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00845 actual_x_data%integral_containers_disk(i)%first%prev => NULL() 00846 actual_x_data%integral_containers_disk(i)%first%next => NULL() 00847 actual_x_data%integral_containers_disk(i)%current => actual_x_data%integral_containers_disk(i)%first 00848 actual_x_data%integral_containers_disk(i)%current%data = 0 00849 actual_x_data%integral_containers_disk(i)%element_counter = 1 00850 actual_x_data%integral_containers_disk(i)%file_counter = 1 00851 actual_x_data%integral_containers_disk(i)%desc = 'Int_' 00852 actual_x_data%integral_containers_disk(i)%unit = -1 00853 WRITE(actual_x_data%integral_containers_disk(i)%filename,'(A,I0,A,A,I0)')& 00854 TRIM(actual_x_data%memory_parameter%storage_location), & 00855 storage_id,"_",actual_x_data%integral_containers_disk(i)%desc, i 00856 CALL compress(actual_x_data%integral_containers_disk(i)%filename, .TRUE.) 00857 END DO 00858 00859 actual_x_data%b_first_load_balance_energy = .TRUE. 00860 actual_x_data%b_first_load_balance_forces = .TRUE. 00861 END DO 00862 END DO 00863 00864 DO irep=1,n_rep_hf 00865 actual_x_data => x_data(irep,1) 00866 CALL hfx_print_info(actual_x_data, hfx_section, irep, error=error) 00867 END DO 00868 CALL timestop(handle) 00869 END SUBROUTINE hfx_create 00870 00871 ! ***************************************************************************** 00880 SUBROUTINE erfc_cutoff(eps,omg,r_cutoff,error) 00881 IMPLICIT NONE 00882 00883 REAL(dp), INTENT(in) :: eps, omg 00884 REAL(dp), INTENT(out) :: r_cutoff 00885 TYPE(cp_error_type), INTENT(inout) :: error 00886 00887 CHARACTER(LEN=*), PARAMETER :: routineN = 'erfc_cutoff', 00888 routineP = moduleN//':'//routineN 00889 00890 REAL(dp), PARAMETER :: abstol=1E-10_dp, soltol=1E-16_dp 00891 REAL(dp) :: r0, f0, fprime0, delta_r 00892 INTEGER :: iter,handle 00893 INTEGER, PARAMETER :: iterMAX = 1000 00894 LOGICAL :: failure 00895 00896 CALL timeset(routineN,handle) 00897 00898 ! initial guess assuming that we are in the asymptotic regime of the erf, and the solution is about 10. 00899 r0 = SQRT(-LOG(eps*omg*10**2))/omg 00900 CALL eval_transc_func(r0,eps,omg,f0,fprime0) 00901 00902 DO iter = 1,iterMAX 00903 delta_r = f0/fprime0 00904 r0 = r0-delta_r 00905 CALL eval_transc_func(r0,eps,omg,f0,fprime0) 00906 IF (ABS(delta_r) .LT. abstol .OR. ABS(f0) .LT. soltol) EXIT 00907 END DO 00908 CPPostcondition(iter<=itermax,cp_failure_level,routineP,error,failure) 00909 r_cutoff = r0 00910 00911 CALL timestop(handle) 00912 CONTAINS 00913 SUBROUTINE eval_transc_func(r,eps,omega,fn,df) 00914 REAL(dp), INTENT(in) :: r, eps, omega 00915 REAL(dp), INTENT(out) :: fn, df 00916 00917 fn = erfc(omega*r) - r*eps 00918 df = -omega*2*EXP(-omega**2*r**2)/SQRT(pi) - eps 00919 END SUBROUTINE eval_transc_func 00920 END SUBROUTINE erfc_cutoff 00921 ! ***************************************************************************** 00926 SUBROUTINE hfx_create_basis_types(basis_parameter,basis_info, atomic_kind_set, do_admm, error) 00927 TYPE(hfx_basis_type), DIMENSION(:), 00928 POINTER :: basis_parameter 00929 TYPE(hfx_basis_info_type) :: basis_info 00930 TYPE(atomic_kind_type), DIMENSION(:), 00931 POINTER :: atomic_kind_set 00932 LOGICAL :: do_admm 00933 TYPE(cp_error_type), INTENT(inout) :: error 00934 00935 CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_create_basis_types', 00936 routineP = moduleN//':'//routineN 00937 00938 INTEGER :: co_counter, handle, i, ikind, ipgf, iset, j, k, la, 00939 max_am_kind, max_coeff, max_nsgfl, max_pgf, max_pgf_kind, max_set, 00940 nkind, nl_count, nset, nseta, offset_a, offset_a1, s_offset_nl_a, sgfa, 00941 so_counter, stat 00942 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nshell 00943 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, nl_a 00944 LOGICAL :: failure 00945 REAL(dp), DIMENSION(:, :), POINTER :: sphi_a 00946 TYPE(atomic_kind_type), POINTER :: atom_kind 00947 TYPE(gto_basis_set_type), POINTER :: orb_basis_a 00948 00949 CALL timeset(routineN,handle) 00950 00951 failure = .FALSE. 00952 00953 !! BASIS parameter 00954 nkind = SIZE(atomic_kind_set,1) 00955 ALLOCATE(basis_parameter(nkind), STAT=stat) 00956 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00957 max_set = 0 00958 DO ikind = 1,nkind 00959 NULLIFY(atom_kind) 00960 atom_kind => atomic_kind_set(ikind) 00961 IF( .NOT. do_admm ) THEN 00962 CALL get_atomic_kind(atomic_kind=atom_kind,& 00963 orb_basis_set=orb_basis_a) 00964 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 00965 maxsgf=basis_info%max_sgf,& 00966 maxnset=basis_info%max_set,& 00967 maxlgto=basis_info%max_am) 00968 00969 ELSE 00970 CALL get_atomic_kind(atomic_kind=atom_kind,& 00971 aux_fit_basis_set=orb_basis_a) 00972 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 00973 maxsgf=basis_info%max_sgf,& 00974 maxnset=basis_info%max_set,& 00975 maxlgto=basis_info%max_am,& 00976 basis_set_id=use_aux_fit_basis_set) 00977 END IF 00978 IF (basis_info%max_set<max_set) STOP "UNEXPECTED MAX_SET" 00979 max_set = MAX(max_set,basis_info%max_set) 00980 CALL get_gto_basis_set(gto_basis_set=orb_basis_a, & 00981 lmax=basis_parameter(ikind)%lmax, & 00982 lmin=basis_parameter(ikind)%lmin, & 00983 npgf=basis_parameter(ikind)%npgf, & 00984 nset=basis_parameter(ikind)%nset, & 00985 zet=basis_parameter(ikind)%zet, & 00986 nsgf_set=basis_parameter(ikind)%nsgf, & 00987 first_sgf=basis_parameter(ikind)%first_sgf, & 00988 sphi=basis_parameter(ikind)%sphi, & 00989 nsgf=basis_parameter(ikind)%nsgf_total,& 00990 l=basis_parameter(ikind)%nl,& 00991 nshell=basis_parameter(ikind)%nshell,& 00992 set_radius=basis_parameter(ikind)%set_radius,& 00993 pgf_radius=basis_parameter(ikind)%pgf_radius,& 00994 kind_radius=basis_parameter(ikind)%kind_radius) 00995 END DO 00996 DO ikind = 1,nkind 00997 ALLOCATE(basis_parameter(ikind)%nsgfl(0:basis_info%max_am,max_set)) 00998 basis_parameter(ikind)%nsgfl = 0 00999 nset = basis_parameter(ikind)%nset 01000 nshell => basis_parameter(ikind)%nshell 01001 DO iset = 1,nset 01002 DO i=0,basis_info%max_am 01003 nl_count = 0 01004 DO j=1,nshell(iset) 01005 IF( basis_parameter(ikind)%nl(j,iset) == i) nl_count = nl_count + 1 01006 END DO 01007 basis_parameter(ikind)%nsgfl(i,iset) = nl_count 01008 END DO 01009 END DO 01010 END DO 01011 01012 max_nsgfl = 0 01013 max_pgf = 0 01014 DO ikind = 1,nkind 01015 max_coeff = 0 01016 max_am_kind = 0 01017 max_pgf_kind = 0 01018 npgfa => basis_parameter(ikind)%npgf 01019 nseta = basis_parameter(ikind)%nset 01020 nl_a => basis_parameter(ikind)%nsgfl 01021 la_max => basis_parameter(ikind)%lmax 01022 la_min => basis_parameter(ikind)%lmin 01023 DO iset = 1,nseta 01024 max_pgf_kind = MAX(max_pgf_kind,npgfa(iset)) 01025 max_pgf = MAX(max_pgf, npgfa(iset)) 01026 DO la = la_min(iset), la_max(iset) 01027 max_nsgfl = MAX(max_nsgfl, nl_a(la,iset)) 01028 max_coeff = MAX(max_coeff, nso(la)*nl_a(la,iset)*nco(la)) 01029 max_am_kind = MAX(max_am_kind, la) 01030 END DO 01031 END DO 01032 ALLOCATE(basis_parameter(ikind)%sphi_ext(max_coeff,0:max_am_kind, max_pgf_kind, nseta), STAT=stat) 01033 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01034 basis_parameter(ikind)%sphi_ext = 0.0_dp 01035 END DO 01036 01037 DO ikind = 1,nkind 01038 sphi_a => basis_parameter(ikind)%sphi 01039 nseta = basis_parameter(ikind)%nset 01040 la_max => basis_parameter(ikind)%lmax 01041 la_min => basis_parameter(ikind)%lmin 01042 npgfa => basis_parameter(ikind)%npgf 01043 first_sgfa => basis_parameter(ikind)%first_sgf 01044 nl_a => basis_parameter(ikind)%nsgfl 01045 DO iset = 1,nseta 01046 sgfa = first_sgfa(1,iset) 01047 DO ipgf = 1,npgfa(iset) 01048 offset_a1 = (ipgf-1)*ncoset(la_max(iset)) 01049 s_offset_nl_a = 0 01050 DO la = la_min(iset),la_max(iset) 01051 offset_a = offset_a1 + ncoset(la-1) 01052 co_counter = 0 01053 co_counter = co_counter + 1 01054 so_counter = 0 01055 DO k = sgfa+s_offset_nl_a,sgfa+s_offset_nl_a+nso(la)*nl_a(la,iset)-1 01056 DO i = offset_a+1,offset_a+nco(la) 01057 so_counter = so_counter + 1 01058 basis_parameter(ikind)%sphi_ext(so_counter,la,ipgf,iset) = sphi_a(i,k) 01059 END DO 01060 END DO 01061 s_offset_nl_a = s_offset_nl_a + nso(la) * (nl_a(la,iset)) 01062 END DO 01063 END DO 01064 END DO 01065 END DO 01066 01067 CALL timestop(handle) 01068 01069 END SUBROUTINE hfx_create_basis_types 01070 01071 01072 SUBROUTINE hfx_release_basis_types(basis_parameter,error) 01073 TYPE(hfx_basis_type), DIMENSION(:), 01074 POINTER :: basis_parameter 01075 TYPE(cp_error_type), INTENT(inout) :: error 01076 01077 CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_release_basis_types', 01078 routineP = moduleN//':'//routineN 01079 01080 INTEGER :: handle, i, stat 01081 LOGICAL :: failure 01082 01083 CALL timeset(routineN,handle) 01084 failure = .FALSE. 01085 01086 !! BASIS parameter 01087 DO i=1,SIZE(basis_parameter) 01088 DEALLOCATE(basis_parameter(i)%nsgfl,STAT=stat) 01089 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01090 DEALLOCATE(basis_parameter(i)%sphi_ext,STAT=stat) 01091 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01092 END DO 01093 DEALLOCATE(basis_parameter,STAT=stat) 01094 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01095 CALL timestop(handle) 01096 01097 END SUBROUTINE hfx_release_basis_types 01098 01099 01100 ! ***************************************************************************** 01103 SUBROUTINE parse_memory_section(memory_parameter, hf_sub_section, storage_id,& 01104 i_thread, n_threads, para_env, irep, skip_disk, skip_in_core_forces, error) 01105 TYPE(hfx_memory_type) :: memory_parameter 01106 TYPE(section_vals_type), POINTER :: hf_sub_section 01107 INTEGER, INTENT(OUT), OPTIONAL :: storage_id 01108 INTEGER, INTENT(IN), OPTIONAL :: i_thread, n_threads 01109 TYPE(cp_para_env_type), OPTIONAL :: para_env 01110 INTEGER, INTENT(IN), OPTIONAL :: irep 01111 LOGICAL, INTENT(IN) :: skip_disk, skip_in_core_forces 01112 TYPE(cp_error_type), INTENT(inout) :: error 01113 01114 CHARACTER(LEN=*), PARAMETER :: routineN = 'parse_memory_section', 01115 routineP = moduleN//':'//routineN 01116 01117 CHARACTER(LEN=512) :: error_msg 01118 CHARACTER(LEN=default_path_length) :: char_val, filename, orig_wd 01119 INTEGER :: int_val, stat 01120 LOGICAL :: check, failure, logic_val 01121 REAL(dp) :: real_val 01122 01123 failure = .FALSE. 01124 check = (PRESENT(storage_id).EQV.PRESENT(i_thread)) .AND.& 01125 (PRESENT(storage_id).EQV.PRESENT(n_threads)).AND.& 01126 (PRESENT(storage_id).EQV.PRESENT(para_env)) .AND.& 01127 (PRESENT(storage_id).EQV.PRESENT(irep)) 01128 CPPostcondition(check,cp_failure_level,routineP,error,failure) 01129 01130 ! Memory Storage 01131 CALL section_vals_val_get(hf_sub_section, "MAX_MEMORY", i_val=int_val, error=error) 01132 memory_parameter%max_memory = int_val 01133 memory_parameter%max_compression_counter = int_val*1024_int_8*128_int_8 01134 CALL section_vals_val_get(hf_sub_section, "EPS_STORAGE", r_val=real_val, error=error) 01135 memory_parameter%eps_storage_scaling = real_val 01136 IF(int_val == 0) THEN 01137 memory_parameter%do_all_on_the_fly = .TRUE. 01138 ELSE 01139 memory_parameter%do_all_on_the_fly = .FALSE. 01140 END IF 01141 memory_parameter%cache_size = CACHE_SIZE 01142 memory_parameter%bits_max_val = BITS_MAX_VAL 01143 memory_parameter%actual_memory_usage = 1 01144 IF( .NOT. skip_in_core_forces) THEN 01145 CALL section_vals_val_get(hf_sub_section, "TREAT_FORCES_IN_CORE", l_val=logic_val, error=error) 01146 memory_parameter%treat_forces_in_core = logic_val 01147 END IF 01148 01149 ! ** IF MAX_MEM == 0 overwrite this flag to false 01150 IF( memory_parameter%do_all_on_the_fly ) memory_parameter%treat_forces_in_core = .FALSE. 01151 01152 ! Disk Storage 01153 IF (.NOT.skip_disk) THEN 01154 memory_parameter%actual_memory_usage_disk = 1 01155 CALL section_vals_val_get(hf_sub_section, "MAX_DISK_SPACE", i_val=int_val, error=error) 01156 memory_parameter%max_compression_counter_disk = int_val*1024_int_8*128_int_8 01157 IF( int_val == 0 ) THEN 01158 memory_parameter%do_disk_storage = .FALSE. 01159 ELSE 01160 memory_parameter%do_disk_storage = .TRUE. 01161 END IF 01162 CALL section_vals_val_get(hf_sub_section, "STORAGE_LOCATION", c_val=char_val, error=error) 01163 CALL compress(char_val,.TRUE.) 01164 !! Add ending / if necessary 01165 01166 IF( SCAN(char_val,"/",.TRUE.) /= LEN_TRIM(char_val)) THEN 01167 WRITE(filename,'(A,A)') TRIM(char_val),"/" 01168 CALL compress(filename) 01169 ELSE 01170 filename = TRIM(char_val) 01171 END IF 01172 CALL compress(filename,.TRUE.) 01173 01174 !! quickly check if we can write on storage_location 01175 CALL m_getcwd(orig_wd) 01176 CALL m_chdir(TRIM(filename),stat) 01177 IF(stat/=0) THEN 01178 WRITE(error_msg,'(A,A,A)') "Request for disk storage failed due to unknown error while writing to ",& 01179 TRIM(filename), ". Please check STORAGE_LOCATION" 01180 CALL cp_assert( .FALSE. , cp_failure_level,cp_assertion_failed,routineP,& 01181 error_msg , & 01182 error,failure) 01183 END IF 01184 CALL m_chdir(orig_wd,stat) 01185 01186 memory_parameter%storage_location = filename 01187 CALL compress(memory_parameter%storage_location, .TRUE.) 01188 ELSE 01189 memory_parameter%do_disk_storage = .FALSE. 01190 END IF 01191 IF (PRESENT(storage_id)) THEN 01192 storage_id = (irep-1)*para_env%num_pe*n_threads + para_env%mepos * n_threads + i_thread-1 01193 END IF 01194 END SUBROUTINE parse_memory_section 01195 01196 ! ***************************************************************************** 01205 SUBROUTINE hfx_release(x_data, error) 01206 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data 01207 TYPE(cp_error_type), INTENT(inout) :: error 01208 01209 CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_release', 01210 routineP = moduleN//':'//routineN 01211 01212 INTEGER :: i, i_thread, irep, n_rep_hf, 01213 n_threads, stat 01214 LOGICAL :: failure 01215 TYPE(hfx_type), POINTER :: actual_x_data 01216 01217 !$ INTEGER :: omp_get_max_threads 01218 01219 01220 failure = .FALSE. 01221 !! There might be 2 hf sections 01222 n_rep_hf = x_data(1,1)%n_rep_hf 01223 n_threads = 1 01224 !$ n_threads = omp_get_max_threads() 01225 01226 IF(x_data(1,1)%potential_parameter%potential_type == do_hfx_potential_truncated .OR. & 01227 x_data(1,1)%potential_parameter%potential_type == do_hfx_potential_mix_cl_trunc ) THEN 01228 init_t_c_g0_lmax = -1 01229 CALL free() 01230 END IF 01231 DO i_thread= 1,n_threads 01232 DO irep=1,n_rep_hf 01233 actual_x_data => x_data(irep,i_thread) 01234 DEALLOCATE(actual_x_data%neighbor_cells,STAT=stat) 01235 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01236 DEALLOCATE(actual_x_data%distribution_energy,STAT=stat) 01237 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01238 DEALLOCATE(actual_x_data%distribution_forces,STAT=stat) 01239 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01240 01241 IF( actual_x_data%load_balance_parameter%blocks_initialized ) THEN 01242 DEALLOCATE(actual_x_data%blocks,STAT=stat) 01243 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01244 IF( i_thread == 1 ) THEN 01245 DEALLOCATE(actual_x_data%pmax_block,STAT=stat) 01246 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01247 END IF 01248 END IF 01249 01250 IF( i_thread == 1 ) THEN 01251 DEALLOCATE(actual_x_data%atomic_pair_list, STAT=stat) 01252 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01253 END IF 01254 01255 IF(actual_x_data%screening_parameter%do_initial_p_screening .OR. & 01256 actual_x_data%screening_parameter%do_p_screening_forces) THEN 01257 IF( i_thread == 1 ) THEN 01258 DEALLOCATE(actual_x_data%pmax_atom, STAT=stat) 01259 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01260 DO i=1,SIZE(actual_x_data%initial_p) 01261 DEALLOCATE(actual_x_data%initial_p(i)%p_kind,STAT=stat) 01262 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01263 END DO 01264 DEALLOCATE(actual_x_data%initial_p,STAT=stat) 01265 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01266 IF( actual_x_data%memory_parameter%treat_forces_in_core) THEN 01267 DEALLOCATE(actual_x_data%pmax_atom_forces, STAT=stat) 01268 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01269 DO i=1,SIZE(actual_x_data%initial_p_forces) 01270 DEALLOCATE(actual_x_data%initial_p_forces(i)%p_kind,STAT=stat) 01271 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01272 END DO 01273 DEALLOCATE(actual_x_data%initial_p_forces,STAT=stat) 01274 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01275 END IF 01276 END IF 01277 DEALLOCATE(actual_x_data%map_atom_to_kind_atom,STAT=stat) 01278 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01279 END IF 01280 IF( i_thread == 1 ) THEN 01281 DEALLOCATE(actual_x_data%is_assoc_atomic_block,STAT=stat) 01282 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01283 DEALLOCATE(actual_x_data%atomic_block_offset,STAT=stat) 01284 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01285 DEALLOCATE(actual_x_data%set_offset,STAT=stat) 01286 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01287 DEALLOCATE(actual_x_data%block_offset,STAT=stat) 01288 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01289 END IF 01290 01291 !! BASIS parameter 01292 CALL hfx_release_basis_types(actual_x_data%basis_parameter,error) 01293 01294 !MK Release libint and libderiv data structure 01295 CALL terminate_libint(actual_x_data%lib) 01296 CALL terminate_libderiv(actual_x_data%lib_deriv) 01297 01298 !! Deallocate containers 01299 CALL dealloc_containers(actual_x_data, hfx_do_eval_energy, error) 01300 CALL dealloc_containers(actual_x_data, hfx_do_eval_forces, error) 01301 01302 !! Deallocate containers 01303 CALL hfx_init_container(actual_x_data%maxval_container_disk, & 01304 actual_x_data%memory_parameter%actual_memory_usage_disk, & 01305 .FALSE., error) 01306 IF( actual_x_data%memory_parameter%do_disk_storage) THEN 01307 CALL close_file(unit_number=actual_x_data%maxval_container_disk%unit,file_status="DELETE") 01308 END IF 01309 DEALLOCATE(actual_x_data%maxval_container_disk%first,STAT=stat) 01310 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01311 DEALLOCATE(actual_x_data%maxval_container_disk,STAT=stat) 01312 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01313 01314 DO i=1,64 01315 CALL hfx_init_container(actual_x_data%integral_containers_disk(i), & 01316 actual_x_data%memory_parameter%actual_memory_usage_disk, & 01317 .FALSE., error) 01318 IF( actual_x_data%memory_parameter%do_disk_storage) THEN 01319 CALL close_file(unit_number=actual_x_data%integral_containers_disk(i)%unit,file_status="DELETE") 01320 END IF 01321 DEALLOCATE(actual_x_data%integral_containers_disk(i)%first,STAT=stat) 01322 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01323 END DO 01324 DEALLOCATE(actual_x_data%integral_containers_disk,STAT=stat) 01325 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01326 01327 01328 ! ** screening functions 01329 IF(actual_x_data%screen_funct_is_initialized) THEN 01330 DEALLOCATE(actual_x_data%screen_funct_coeffs_set,STAT=stat) 01331 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01332 DEALLOCATE(actual_x_data%screen_funct_coeffs_kind,STAT=stat) 01333 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01334 DEALLOCATE(actual_x_data%pair_dist_radii_pgf,STAT=stat) 01335 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01336 DEALLOCATE(actual_x_data%screen_funct_coeffs_pgf,STAT=stat) 01337 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01338 actual_x_data%screen_funct_is_initialized = .FALSE. 01339 END IF 01340 01341 ! ** maps 01342 IF( ASSOCIATED(actual_x_data%map_atoms_to_cpus) ) THEN 01343 DO i=1,SIZE(actual_x_data%map_atoms_to_cpus) 01344 DEALLOCATE(actual_x_data%map_atoms_to_cpus(i)%iatom_list,STAT=stat) 01345 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01346 DEALLOCATE(actual_x_data%map_atoms_to_cpus(i)%jatom_list,STAT=stat) 01347 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01348 END DO 01349 DEALLOCATE(actual_x_data%map_atoms_to_cpus,STAT=stat) 01350 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01351 END IF 01352 END DO 01353 01354 END DO 01355 DEALLOCATE(x_data,STAT=stat) 01356 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01357 END SUBROUTINE hfx_release 01358 01359 ! ***************************************************************************** 01373 SUBROUTINE hfx_create_neighbor_cells(x_data, pbc_shells, cell, i_thread, & 01374 error) 01375 TYPE(hfx_type), POINTER :: x_data 01376 INTEGER, INTENT(INOUT) :: pbc_shells 01377 TYPE(cell_type), POINTER :: cell 01378 INTEGER, INTENT(IN) :: i_thread 01379 TYPE(cp_error_type), INTENT(inout) :: error 01380 01381 CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_create_neighbor_cells', 01382 routineP = moduleN//':'//routineN 01383 01384 CHARACTER(LEN=512) :: error_msg 01385 CHARACTER(LEN=64) :: char_nshells 01386 INTEGER :: i, idx, ikind, ipgf, iset, ishell, j, jkind, jpgf, jset, 01387 jshell, k, kshell, l, m(3), max_shell, nseta, nsetb, stat, 01388 total_number_of_cells, ub, ub_max 01389 INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, npgfa, npgfb 01390 LOGICAL :: failure, image_cell_found, 01391 nothing_more_to_add 01392 REAL(dp) :: cross_product(3), dist_min, distance(14), l_min, normal(3,6), 01393 P(3,14), plane_vector(3,2), point_in_plane(3), r(3), R1, R_max, 01394 R_max_stress, s(3), x, y, z, Zeta1 01395 REAL(dp), DIMENSION(:, :), POINTER :: zeta, zetb 01396 TYPE(hfx_cell_type), ALLOCATABLE, 01397 DIMENSION(:) :: tmp_neighbor_cells 01398 01399 failure = .FALSE. 01400 total_number_of_cells = 0 01401 01402 ! ** Check some settings 01403 IF( i_thread == 1 ) THEN 01404 IF( x_data%potential_parameter%potential_type /= do_hfx_potential_truncated .AND. & 01405 x_data%potential_parameter%potential_type /= do_hfx_potential_short .AND. & 01406 x_data%potential_parameter%potential_type /= do_hfx_potential_mix_cl_trunc) THEN 01407 CALL cp_assert(.FALSE.,cp_warning_level,cp_assertion_failed,routineP,& 01408 "Periodic Hartree Fock calculation requested without use "//& 01409 "of a truncated or shortrange potential. This may lead to unphysical total energies. "//& 01410 "Use a truncated potential to avoid possible problems. "//& 01411 CPSourceFileRef,& 01412 only_ionode=.TRUE.) 01413 ELSE 01414 l_min = MIN( plane_distance(1,0,0,cell), & 01415 plane_distance(0,1,0,cell), & 01416 plane_distance(0,0,1,cell) ) 01417 l_min = 0.5_dp * l_min 01418 IF( x_data%potential_parameter%cutoff_radius >= l_min ) THEN 01419 CALL cp_assert(.FALSE.,cp_warning_level,cp_assertion_failed,routineP,& 01420 "Periodic Hartree Fock calculation requested with use "//& 01421 "of a truncated or shortrange potential. The cutoff radius is larger than half "//& 01422 "the minimal cell dimension. This may lead to unphysical "//& 01423 "total energies. Reduce the cutoff radius in order to avoid "//& 01424 "possible problems. "//& 01425 CPSourceFileRef,& 01426 only_ionode=.TRUE.) 01427 END IF 01428 END IF 01429 END IF 01430 01431 SELECT CASE(x_data%potential_parameter%potential_type) 01432 CASE( do_hfx_potential_truncated, do_hfx_potential_mix_cl_trunc, do_hfx_potential_short ) 01433 R_max = 0.0_dp 01434 DO ikind=1,SIZE(x_data%basis_parameter) 01435 la_max => x_data%basis_parameter(ikind)%lmax 01436 zeta => x_data%basis_parameter(ikind)%zet 01437 nseta = x_data%basis_parameter(ikind)%nset 01438 npgfa => x_data%basis_parameter(ikind)%npgf 01439 DO jkind=1,SIZE(x_data%basis_parameter) 01440 lb_max => x_data%basis_parameter(jkind)%lmax 01441 zetb => x_data%basis_parameter(jkind)%zet 01442 nsetb = x_data%basis_parameter(jkind)%nset 01443 npgfb => x_data%basis_parameter(jkind)%npgf 01444 DO iset=1,nseta 01445 DO jset=1,nsetb 01446 DO ipgf=1,npgfa(iset) 01447 DO jpgf=1,npgfb(jset) 01448 Zeta1 = zeta(ipgf,iset) + zetb(jpgf,jset) 01449 R1 = 1.0_dp/SQRT(Zeta1) * mul_fact(la_max(iset)+lb_max(jset)) * & 01450 SQRT(-LOG(x_data%screening_parameter%eps_schwarz)) 01451 R_max = MAX(R1,R_max) 01452 END DO 01453 END DO 01454 END DO 01455 END DO 01456 END DO 01457 END DO 01458 01459 R_max = 2.0_dp*R_max + x_data%potential_parameter%cutoff_radius 01460 nothing_more_to_add = .FALSE. 01461 max_shell = 0 01462 total_number_of_cells = 0 01463 ub = 1 01464 DEALLOCATE(x_data%neighbor_cells,STAT=stat) 01465 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01466 ALLOCATE(x_data%neighbor_cells(1), STAT=stat) 01467 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01468 x_data%neighbor_cells(1)%cell = 0.0_dp 01469 x_data%neighbor_cells(1)%cell_r = 0.0_dp 01470 01471 ! ** What follows is kind of a ray tracing algorithm 01472 ! ** Given a image cell (ishell, jshell, kshell) we try to figure out the 01473 ! ** shortest distance of this image cell to the basic unit cell (0,0,0), i.e. the point 01474 ! ** (0.0, 0.0, 0.0) 01475 ! ** This is achieved by checking the 8 Corners of the cell, and, in addition, the shortest distance 01476 ! ** to all 6 faces. The faces are only taken into account if the penetration point of the normal 01477 ! ** to the plane defined by a face lies within this face. 01478 ! ** This is very fast, because no trigonometric functions are being used 01479 ! ** The points are defined as follows 01480 ! ** 01481 ! ** 01482 ! ** _________________________ 01483 ! ** /P4____________________P8/| 01484 ! ** / / ___________________/ / | 01485 ! ** / / /| | / / | z 01486 ! ** / / / | | / / . | /|\ _ y 01487 ! ** / / /| | | / / /| | | /| 01488 ! ** / / / | | | / / / | | | / 01489 ! ** / / / | | | / / /| | | | / 01490 ! ** / /_/___| | |__________/ / / | | | |/ 01491 ! ** /P2______| | |_________P6/ / | | | ----------> x 01492 ! ** | _______| | |_________| | | | | | 01493 ! ** | | | | | |________________| | | 01494 ! ** | | | |P3___________________P7 | 01495 ! ** | | | / / _________________ / / 01496 ! ** | | | / / / | | |/ / / 01497 ! ** | | | / / / | | | / / 01498 ! ** | | |/ / / | | |/ / 01499 ! ** | | | / / | | ' / 01500 ! ** | | |/_/_______________| | / 01501 ! ** | |____________________| | / 01502 ! ** |P1_____________________P5/ 01503 ! ** 01504 ! ** 01505 01506 DO WHILE (.NOT. nothing_more_to_add) 01507 ! Calculate distances to the eight points P1 to P8 01508 image_cell_found = .FALSE. 01509 ALLOCATE(tmp_neighbor_cells(1:ub),STAT=stat) 01510 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01511 DO i = 1,ub-1 01512 tmp_neighbor_cells(i) = x_data%neighbor_cells(i) 01513 END DO 01514 ub_max = (2*max_shell+1)**3 01515 DEALLOCATE(x_data%neighbor_cells, STAT=stat) 01516 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01517 ALLOCATE(x_data%neighbor_cells(1:ub_max), STAT=stat) 01518 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01519 DO i = 1,ub-1 01520 x_data%neighbor_cells(i) = tmp_neighbor_cells(i) 01521 END DO 01522 DO i = ub,ub_max 01523 x_data%neighbor_cells(i)%cell = 0.0_dp 01524 x_data%neighbor_cells(i)%cell_r = 0.0_dp 01525 END DO 01526 01527 DEALLOCATE(tmp_neighbor_cells, STAT=stat) 01528 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01529 01530 DO ishell=-max_shell,max_shell 01531 DO jshell=-max_shell,max_shell 01532 DO kshell=-max_shell,max_shell 01533 IF( MAX(ABS(ishell),ABS(jshell),ABS(kshell)) /= max_shell) CYCLE 01534 idx = 0 01535 DO j=0,1 01536 x = -1.0_dp/2.0_dp + j*1.0_dp 01537 DO k=0,1 01538 y = -1.0_dp/2.0_dp + k*1.0_dp 01539 DO l=0,1 01540 z = -1.0_dp/2.0_dp + l*1.0_dp 01541 idx = idx + 1 01542 P(1,idx) = x + ishell 01543 P(2,idx) = y + jshell 01544 P(3,idx) = z + kshell 01545 CALL scaled_to_real(r,P(:,idx),cell) 01546 distance(idx) = SQRT(SUM(r**2)) 01547 P(1:3,idx)=r 01548 END DO 01549 END DO 01550 END DO 01551 ! Now check distance to Faces and only take them into account if the base point lies within quadrilateral 01552 01553 ! Face A (1342) 1 is the reference 01554 idx = idx + 1 01555 plane_vector(:,1) = P(:,3)-P(:,1) 01556 plane_vector(:,2) = P(:,2)-P(:,1) 01557 cross_product(1) = plane_vector(2,1)*plane_vector(3,2) - plane_vector(3,1)*plane_vector(2,2) 01558 cross_product(2) = plane_vector(3,1)*plane_vector(1,2) - plane_vector(1,1)*plane_vector(3,2) 01559 cross_product(3) = plane_vector(1,1)*plane_vector(2,2) - plane_vector(2,1)*plane_vector(1,2) 01560 normal(:,1) = cross_product/SQRT(SUM(cross_product**2)) 01561 point_in_plane = -normal(:,1)*(normal(1,1)*P(1,1)+normal(2,1)*P(2,1)+normal(3,1)*P(3,1)) 01562 01563 IF (point_is_in_quadrilateral(P(:,1),P(:,3),P(:,4),P(:,2),point_in_plane) ) THEN 01564 distance(idx)=ABS(normal(1,1)*P(1,1)+normal(2,1)*P(2,1)+normal(3,1)*P(3,1)) 01565 ELSE 01566 distance(idx)=HUGE(distance(idx)) 01567 END IF 01568 01569 ! Face B (1562) 1 is the reference 01570 idx = idx + 1 01571 plane_vector(:,1) = P(:,2)-P(:,1) 01572 plane_vector(:,2) = P(:,5)-P(:,1) 01573 cross_product(1) = plane_vector(2,1)*plane_vector(3,2) - plane_vector(3,1)*plane_vector(2,2) 01574 cross_product(2) = plane_vector(3,1)*plane_vector(1,2) - plane_vector(1,1)*plane_vector(3,2) 01575 cross_product(3) = plane_vector(1,1)*plane_vector(2,2) - plane_vector(2,1)*plane_vector(1,2) 01576 normal(:,1) = cross_product/SQRT(SUM(cross_product**2)) 01577 point_in_plane = -normal(:,1)*(normal(1,1)*P(1,1)+normal(2,1)*P(2,1)+normal(3,1)*P(3,1)) 01578 01579 IF (point_is_in_quadrilateral(P(:,1),P(:,5),P(:,6),P(:,2),point_in_plane) ) THEN 01580 distance(idx)=ABS(normal(1,1)*P(1,1)+normal(2,1)*P(2,1)+normal(3,1)*P(3,1)) 01581 ELSE 01582 distance(idx)=HUGE(distance(idx)) 01583 END IF 01584 01585 ! Face C (5786) 5 is the reference 01586 idx = idx + 1 01587 plane_vector(:,1) = P(:,7)-P(:,5) 01588 plane_vector(:,2) = P(:,6)-P(:,5) 01589 cross_product(1) = plane_vector(2,1)*plane_vector(3,2) - plane_vector(3,1)*plane_vector(2,2) 01590 cross_product(2) = plane_vector(3,1)*plane_vector(1,2) - plane_vector(1,1)*plane_vector(3,2) 01591 cross_product(3) = plane_vector(1,1)*plane_vector(2,2) - plane_vector(2,1)*plane_vector(1,2) 01592 normal(:,1) = cross_product/SQRT(SUM(cross_product**2)) 01593 point_in_plane = -normal(:,1)*(normal(1,1)*P(1,5)+normal(2,1)*P(2,5)+normal(3,1)*P(3,5)) 01594 01595 IF (point_is_in_quadrilateral(P(:,5),P(:,7),P(:,8),P(:,6),point_in_plane) ) THEN 01596 distance(idx)=ABS(normal(1,1)*P(1,5)+normal(2,1)*P(2,5)+normal(3,1)*P(3,5)) 01597 ELSE 01598 distance(idx)=HUGE(distance(idx)) 01599 END IF 01600 01601 ! Face D (3784) 3 is the reference 01602 idx = idx + 1 01603 plane_vector(:,1) = P(:,7)-P(:,3) 01604 plane_vector(:,2) = P(:,4)-P(:,3) 01605 cross_product(1) = plane_vector(2,1)*plane_vector(3,2) - plane_vector(3,1)*plane_vector(2,2) 01606 cross_product(2) = plane_vector(3,1)*plane_vector(1,2) - plane_vector(1,1)*plane_vector(3,2) 01607 cross_product(3) = plane_vector(1,1)*plane_vector(2,2) - plane_vector(2,1)*plane_vector(1,2) 01608 normal(:,1) = cross_product/SQRT(SUM(cross_product**2)) 01609 point_in_plane = -normal(:,1)*(normal(1,1)*P(1,3)+normal(2,1)*P(2,3)+normal(3,1)*P(3,3)) 01610 01611 IF (point_is_in_quadrilateral(P(:,3),P(:,7),P(:,8),P(:,4),point_in_plane) ) THEN 01612 distance(idx)=ABS(normal(1,1)*P(1,3)+normal(2,1)*P(2,3)+normal(3,1)*P(3,3)) 01613 ELSE 01614 distance(idx)=HUGE(distance(idx)) 01615 END IF 01616 01617 ! Face E (2684) 2 is the reference 01618 idx = idx + 1 01619 plane_vector(:,1) = P(:,6)-P(:,2) 01620 plane_vector(:,2) = P(:,4)-P(:,2) 01621 cross_product(1) = plane_vector(2,1)*plane_vector(3,2) - plane_vector(3,1)*plane_vector(2,2) 01622 cross_product(2) = plane_vector(3,1)*plane_vector(1,2) - plane_vector(1,1)*plane_vector(3,2) 01623 cross_product(3) = plane_vector(1,1)*plane_vector(2,2) - plane_vector(2,1)*plane_vector(1,2) 01624 normal(:,1) = cross_product/SQRT(SUM(cross_product**2)) 01625 point_in_plane = -normal(:,1)*(normal(1,1)*P(1,2)+normal(2,1)*P(2,2)+normal(3,1)*P(3,2)) 01626 01627 IF (point_is_in_quadrilateral(P(:,2),P(:,6),P(:,8),P(:,4),point_in_plane) ) THEN 01628 distance(idx)=ABS(normal(1,1)*P(1,2)+normal(2,1)*P(2,2)+normal(3,1)*P(3,2)) 01629 ELSE 01630 distance(idx)=HUGE(distance(idx)) 01631 END IF 01632 01633 ! Face F (1573) 1 is the reference 01634 idx = idx + 1 01635 plane_vector(:,1) = P(:,5)-P(:,1) 01636 plane_vector(:,2) = P(:,3)-P(:,1) 01637 cross_product(1) = plane_vector(2,1)*plane_vector(3,2) - plane_vector(3,1)*plane_vector(2,2) 01638 cross_product(2) = plane_vector(3,1)*plane_vector(1,2) - plane_vector(1,1)*plane_vector(3,2) 01639 cross_product(3) = plane_vector(1,1)*plane_vector(2,2) - plane_vector(2,1)*plane_vector(1,2) 01640 normal(:,1) = cross_product/SQRT(SUM(cross_product**2)) 01641 point_in_plane = -normal(:,1)*(normal(1,1)*P(1,1)+normal(2,1)*P(2,1)+normal(3,1)*P(3,1)) 01642 01643 IF (point_is_in_quadrilateral(P(:,1),P(:,5),P(:,7),P(:,3),point_in_plane) )THEN 01644 distance(idx)=ABS(normal(1,1)*P(1,1)+normal(2,1)*P(2,1)+normal(3,1)*P(3,1)) 01645 ELSE 01646 distance(idx)=HUGE(distance(idx)) 01647 END IF 01648 01649 dist_min=MINVAL(distance) 01650 IF( max_shell == 0) THEN 01651 image_cell_found=.TRUE. 01652 END IF 01653 IF(dist_min<R_max)THEN 01654 total_number_of_cells = total_number_of_cells + 1 01655 x_data%neighbor_cells(ub)%cell = REAL((/ishell,jshell,kshell/),dp) 01656 ub = ub + 1 01657 image_cell_found = .TRUE. 01658 END IF 01659 01660 END DO 01661 END DO 01662 END DO 01663 IF( image_cell_found) THEN 01664 max_shell = max_shell + 1 01665 ELSE 01666 nothing_more_to_add = .TRUE. 01667 END IF 01668 END DO 01669 ! now remove what is not needed 01670 ALLOCATE(tmp_neighbor_cells(total_number_of_cells), STAT=stat) 01671 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01672 DO i=1,ub-1 01673 tmp_neighbor_cells(i) = x_data%neighbor_cells(i) 01674 END DO 01675 DEALLOCATE(x_data%neighbor_cells, STAT=stat) 01676 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01677 ! If we only need the supercell, total_number_of_cells is still 0, repair 01678 IF( total_number_of_cells == 0 ) THEN 01679 total_number_of_cells = 1 01680 ALLOCATE(x_data%neighbor_cells(total_number_of_cells), STAT=stat) 01681 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01682 DO i=1,total_number_of_cells 01683 x_data%neighbor_cells(i)%cell = 0.0_dp 01684 x_data%neighbor_cells(i)%cell_r = 0.0_dp 01685 END DO 01686 ELSE 01687 ALLOCATE(x_data%neighbor_cells(total_number_of_cells), STAT=stat) 01688 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01689 DO i=1,total_number_of_cells 01690 x_data%neighbor_cells(i) = tmp_neighbor_cells(i) 01691 END DO 01692 END IF 01693 DEALLOCATE(tmp_neighbor_cells, STAT=stat) 01694 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01695 01696 IF(x_data%periodic_parameter%number_of_shells == do_hfx_auto_shells ) THEN 01697 ! Do nothing 01698 ELSE 01699 total_number_of_cells = 0 01700 DO i = 0,x_data%periodic_parameter%number_of_shells 01701 total_number_of_cells = total_number_of_cells + count_cells_perd(i,x_data%periodic_parameter%perd) 01702 END DO 01703 IF( total_number_of_cells < SIZE(x_data%neighbor_cells)) THEN 01704 IF( i_thread == 1) THEN 01705 01706 WRITE(char_nshells,'(I3)') SIZE(x_data%neighbor_cells) 01707 WRITE(error_msg,'(A,A,A)') "Periodic Hartree Fock calculation requested with use "//& 01708 "of a truncated potential. The number of shells to be considered "//& 01709 "might be too small. CP2K conservatively estimates to need "//TRIM(char_nshells)//" periodic images "//& 01710 "Please carefully check if you get converged results." 01711 01712 CALL cp_assert(.FALSE.,cp_warning_level,cp_assertion_failed,routineP,& 01713 error_msg,& 01714 error,& 01715 only_ionode=.TRUE.) 01716 END IF 01717 END IF 01718 total_number_of_cells = 0 01719 DO i = 0,x_data%periodic_parameter%number_of_shells 01720 total_number_of_cells = total_number_of_cells + count_cells_perd(i,x_data%periodic_parameter%perd) 01721 END DO 01722 DEALLOCATE(x_data%neighbor_cells, STAT=stat) 01723 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01724 01725 ALLOCATE(x_data%neighbor_cells(total_number_of_cells), STAT=stat) 01726 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01727 m = 0 01728 i = 1 01729 DO WHILE(SUM(m**2)<=x_data%periodic_parameter%number_of_shells) 01730 x_data%neighbor_cells(i)%cell = REAL(m,dp) 01731 CALL next_image_cell_perd(m,x_data%periodic_parameter%perd) 01732 i=i+1 01733 ENDDO 01734 END IF 01735 CASE DEFAULT 01736 total_number_of_cells = 0 01737 IF( pbc_shells == -1 ) pbc_shells = 0 01738 DO i = 0,pbc_shells 01739 total_number_of_cells = total_number_of_cells + count_cells_perd(i,x_data%periodic_parameter%perd) 01740 END DO 01741 DEALLOCATE(x_data%neighbor_cells,STAT=stat) 01742 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01743 01744 ALLOCATE(x_data%neighbor_cells(total_number_of_cells), STAT=stat) 01745 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01746 01747 m = 0 01748 i = 1 01749 DO WHILE(SUM(m**2)<=pbc_shells) 01750 x_data%neighbor_cells(i)%cell = REAL(m,dp) 01751 CALL next_image_cell_perd(m,x_data%periodic_parameter%perd) 01752 i=i+1 01753 ENDDO 01754 END SELECT 01755 01756 ! ** Transform into real coord 01757 DO i=1,SIZE(x_data%neighbor_cells) 01758 r = 0.0_dp 01759 x_data%neighbor_cells(i)%cell_r(:) = 0.0_dp 01760 s = x_data%neighbor_cells(i)%cell(:) 01761 CALL scaled_to_real(x_data%neighbor_cells(i)%cell_r,s,cell) 01762 END DO 01763 x_data%periodic_parameter%number_of_shells = pbc_shells 01764 01765 R_max_stress = 0.0_dp 01766 DO i=1,SIZE(x_data%neighbor_cells) 01767 R_max_stress = MAX(R_max_stress, MAXVAL(ABS(x_data%neighbor_cells(i)%cell_r(:)))) 01768 END DO 01769 R_max_stress = R_max_stress + ABS(MAXVAL(cell%hmat(:,:))) 01770 x_data%periodic_parameter%R_max_stress = R_max_stress 01771 01772 END SUBROUTINE hfx_create_neighbor_cells 01773 01774 FUNCTION point_is_in_quadrilateral(A,B,C,D,P) 01775 REAL(dp) :: A(3), B(3), C(3), D(3), P(3) 01776 LOGICAL :: point_is_in_quadrilateral 01777 01778 REAL(dp) :: dot00, dot01, dot02, dot11, 01779 dot12, invDenom, u, v, v0(3), 01780 v1(3), v2(3) 01781 01782 point_is_in_quadrilateral = .FALSE. 01783 01784 01785 ! ** Check for both triangles ABC and ACD 01786 ! ** 01787 ! ** D -------------- C 01788 ! ** / / 01789 ! ** / / 01790 ! ** A----------------B 01791 ! ** 01792 ! ** 01793 ! ** 01794 01795 ! ** ABC 01796 01797 v0 = D - A 01798 v1 = C - A 01799 v2 = P - A 01800 01801 ! ** Compute dot products 01802 dot00 = DOT_PRODUCT(v0, v0) 01803 dot01 = DOT_PRODUCT(v0, v1) 01804 dot02 = DOT_PRODUCT(v0, v2) 01805 dot11 = DOT_PRODUCT(v1, v1) 01806 dot12 = DOT_PRODUCT(v1, v2) 01807 01808 ! ** Compute barycentric coordinates 01809 invDenom = 1 / (dot00 * dot11 - dot01 * dot01) 01810 u = (dot11 * dot02 - dot01 * dot12) * invDenom 01811 v = (dot00 * dot12 - dot01 * dot02) * invDenom 01812 ! ** Check if point is in triangle 01813 IF ( (u >= 0) .AND. (v >= 0) .AND. (u + v <= 1) ) THEN 01814 point_is_in_quadrilateral = .TRUE. 01815 RETURN 01816 END IF 01817 v0 = C - A 01818 v1 = B - A 01819 v2 = P - A 01820 01821 ! ** Compute dot products 01822 dot00 = DOT_PRODUCT(v0, v0) 01823 dot01 = DOT_PRODUCT(v0, v1) 01824 dot02 = DOT_PRODUCT(v0, v2) 01825 dot11 = DOT_PRODUCT(v1, v1) 01826 dot12 = DOT_PRODUCT(v1, v2) 01827 01828 ! ** Compute barycentric coordinates 01829 invDenom = 1 / (dot00 * dot11 - dot01 * dot01) 01830 u = (dot11 * dot02 - dot01 * dot12) * invDenom 01831 v = (dot00 * dot12 - dot01 * dot02) * invDenom 01832 01833 ! ** Check if point is in triangle 01834 IF ( (u >= 0) .AND. (v >= 0) .AND. (u + v <= 1) ) THEN 01835 point_is_in_quadrilateral = .TRUE. 01836 RETURN 01837 END IF 01838 01839 END FUNCTION point_is_in_quadrilateral 01840 01841 01842 ! ***************************************************************************** 01852 SUBROUTINE hfx_init_container(container, memory_usage, do_disk_storage, error) 01853 TYPE(hfx_container_type) :: container 01854 INTEGER :: memory_usage 01855 LOGICAL :: do_disk_storage 01856 TYPE(cp_error_type), INTENT(inout) :: error 01857 01858 CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_init_container', 01859 routineP = moduleN//':'//routineN 01860 01861 INTEGER :: stat 01862 LOGICAL :: failure 01863 TYPE(hfx_container_node), POINTER :: current, next 01864 01865 failure =.FALSE. 01866 01867 !! DEALLOCATE memory 01868 01869 current => container%first 01870 DO WHILE(ASSOCIATED(current)) 01871 next => current%next 01872 DEALLOCATE(current) 01873 current=>next 01874 END DO 01875 01876 !! Allocate first list entry, init members 01877 ALLOCATE(container%first,STAT=stat) 01878 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01879 container%first%prev => NULL() 01880 container%first%next => NULL() 01881 container%current => container%first 01882 container%current%data = 0 01883 container%element_counter = 1 01884 memory_usage = 1 01885 01886 IF( do_disk_storage) THEN 01887 !! close the file, if this is no the first time 01888 IF(container%unit /= -1 ) THEN 01889 CALL close_file(unit_number=container%unit) 01890 END IF 01891 CALL open_file(file_name=TRIM(container%filename),file_status="UNKNOWN",file_form="UNFORMATTED",file_action="WRITE",& 01892 unit_number=container%unit) 01893 END IF 01894 01895 END SUBROUTINE hfx_init_container 01896 01897 01898 ! ***************************************************************************** 01909 SUBROUTINE hfx_set_distr_energy(ptr_to_distr,x_data, error) 01910 TYPE(hfx_distribution), DIMENSION(:), 01911 POINTER :: ptr_to_distr 01912 TYPE(hfx_type), POINTER :: x_data 01913 TYPE(cp_error_type), INTENT(inout) :: error 01914 01915 CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_set_distr_energy', 01916 routineP = moduleN//':'//routineN 01917 01918 INTEGER :: stat 01919 LOGICAL :: failure 01920 01921 DEALLOCATE(x_data%distribution_energy,STAT=stat) 01922 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01923 01924 ALLOCATE(x_data%distribution_energy(SIZE(ptr_to_distr)),STAT=stat) 01925 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01926 x_data%distribution_energy = ptr_to_distr 01927 01928 END SUBROUTINE hfx_set_distr_energy 01929 01930 ! ***************************************************************************** 01941 SUBROUTINE hfx_set_distr_forces(ptr_to_distr, x_data, error) 01942 TYPE(hfx_distribution), DIMENSION(:), 01943 POINTER :: ptr_to_distr 01944 TYPE(hfx_type), POINTER :: x_data 01945 TYPE(cp_error_type), INTENT(inout) :: error 01946 01947 CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_set_distr_forces', 01948 routineP = moduleN//':'//routineN 01949 01950 INTEGER :: stat 01951 LOGICAL :: failure 01952 01953 DEALLOCATE(x_data%distribution_forces,STAT=stat) 01954 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01955 01956 ALLOCATE(x_data%distribution_forces(SIZE(ptr_to_distr)),STAT=stat) 01957 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01958 x_data%distribution_forces = ptr_to_distr 01959 01960 END SUBROUTINE hfx_set_distr_forces 01961 01962 ! ***************************************************************************** 01974 SUBROUTINE hfx_reset_memory_usage_counter(memory_parameter, subtr_size_mb) 01975 01976 TYPE(hfx_memory_type) :: memory_parameter 01977 INTEGER(int_8), INTENT(IN) :: subtr_size_mb 01978 01979 INTEGER(int_8) :: max_memory 01980 01981 max_memory = memory_parameter%max_memory 01982 max_memory = max_memory - subtr_size_mb 01983 IF( max_memory <= 0 ) THEN 01984 memory_parameter%do_all_on_the_fly = .TRUE. 01985 memory_parameter%max_compression_counter = 0 01986 ELSE 01987 memory_parameter%do_all_on_the_fly = .FALSE. 01988 memory_parameter%max_compression_counter = max_memory*1024_int_8*128_int_8 01989 END IF 01990 END SUBROUTINE hfx_reset_memory_usage_counter 01991 01992 ! ***************************************************************************** 02003 SUBROUTINE hfx_print_info(x_data, hfx_section, i_rep, error) 02004 TYPE(hfx_type), POINTER :: x_data 02005 TYPE(section_vals_type), POINTER :: hfx_section 02006 INTEGER, INTENT(IN) :: i_rep 02007 TYPE(cp_error_type), INTENT(inout) :: error 02008 02009 INTEGER :: iw 02010 REAL(dp) :: rc_ang 02011 TYPE(cp_logger_type), POINTER :: logger 02012 02013 NULLIFY(logger) 02014 logger => cp_error_get_logger(error) 02015 02016 iw = cp_print_key_unit_nr(logger,hfx_section,"HF_INFO",& 02017 extension=".scfLog",error=error) 02018 IF (iw>0) THEN 02019 WRITE (UNIT=iw,FMT="(/,(T3,A,T61,I20))")& 02020 "HFX_INFO| Replica ID: ",i_rep 02021 WRITE (UNIT=iw,FMT="((T3,A,T73,ES8.1))")& 02022 "HFX_INFO| EPS_SCHWARZ: ", x_data%screening_parameter%eps_schwarz 02023 WRITE (UNIT=iw,FMT="((T3,A,T73,ES8.1))")& 02024 "HFX_INFO| EPS_SCHWARZ_FORCES ", x_data%screening_parameter%eps_schwarz_forces 02025 WRITE (UNIT=iw,FMT="((T3,A,T73,ES8.1))")& 02026 "HFX_INFO| EPS_STORAGE_SCALING: ", x_data%memory_parameter%eps_storage_scaling 02027 WRITE (UNIT=iw,FMT="((T3,A,T61,I20))")& 02028 "HFX_INFO| NBINS: ", x_data%load_balance_parameter%nbins 02029 WRITE (UNIT=iw,FMT="((T3,A,T61,I20))")& 02030 "HFX_INFO| BLOCK_SIZE: ", x_data%load_balance_parameter%block_size 02031 WRITE(iw,'(T3,A,T61,F20.10)') & 02032 "HFX_INFO| FRACTION: ", x_data%general_parameter%fraction 02033 SELECT CASE (x_data%potential_parameter%potential_type) 02034 CASE(do_hfx_potential_coulomb) 02035 WRITE (UNIT=iw,FMT="((T3,A,T74,A))")& 02036 "HFX_INFO| Interaction Potential: ", "COULOMB" 02037 CASE(do_hfx_potential_short) 02038 WRITE (UNIT=iw,FMT="((T3,A,T71,A))")& 02039 "HFX_INFO| Interaction Potential: ", "SHORTRANGE" 02040 WRITE(iw,'(T3,A,T61,F20.10)') & 02041 "HFX_INFO| Omega: ", x_data%potential_parameter%omega 02042 rc_ang = cp_unit_from_cp2k(x_data%potential_parameter%cutoff_radius,"angstrom",error=error) 02043 WRITE(iw,'(T3,A,T61,F20.10)') & 02044 "HFX_INFO| Cutoff Radius [angstrom]: ", rc_ang 02045 CASE(do_hfx_potential_long) 02046 WRITE (UNIT=iw,FMT="((T3,A,T72,A))")& 02047 "HFX_INFO| Interaction Potential: ", "LONGRANGE" 02048 WRITE(iw,'(T3,A,T61,F20.10)') & 02049 "HFX_INFO| Omega: ", x_data%potential_parameter%omega 02050 CASE(do_hfx_potential_mix_cl) 02051 WRITE (UNIT=iw,FMT="((T3,A,T75,A))")& 02052 "HFX_INFO| Interaction Potential: ", "MIX_CL" 02053 WRITE(iw,'(T3,A,T61,F20.10)') & 02054 "HFX_INFO| Omega: ", x_data%potential_parameter%omega 02055 WRITE(iw,'(T3,A,T61,F20.10)') & 02056 "HFX_INFO| SCALE_COULOMB: ", x_data%potential_parameter%scale_coulomb 02057 WRITE(iw,'(T3,A,T61,F20.10)') & 02058 "HFX_INFO| SCALE_LONGRANGE: ", x_data%potential_parameter%scale_longrange 02059 CASE(do_hfx_potential_gaussian) 02060 WRITE (UNIT=iw,FMT="((T3,A,T73,A))")& 02061 "HFX_INFO| Interaction Potential: ", "GAUSSIAN" 02062 WRITE(iw,'(T3,A,T61,F20.10)') & 02063 "HFX_INFO| Omega: ", x_data%potential_parameter%omega 02064 CASE(do_hfx_potential_mix_lg) 02065 WRITE (UNIT=iw,FMT="((T3,A,T75,A))")& 02066 "HFX_INFO| Interaction Potential: ", "MIX_LG" 02067 WRITE(iw,'(T3,A,T61,F20.10)') & 02068 "HFX_INFO| Omega: ", x_data%potential_parameter%omega 02069 WRITE(iw,'(T3,A,T61,F20.10)') & 02070 "HFX_INFO| SCALE_LONGRANGE: ", x_data%potential_parameter%scale_longrange 02071 WRITE(iw,'(T3,A,T61,F20.10)') & 02072 "HFX_INFO| SCALE_GAUSSIAN: ", x_data%potential_parameter%scale_gaussian 02073 CASE(do_hfx_potential_id) 02074 WRITE (UNIT=iw,FMT="((T3,A,T73,A))")& 02075 "HFX_INFO| Interaction Potential: ", "IDENTITY" 02076 CASE(do_hfx_potential_truncated) 02077 WRITE (UNIT=iw,FMT="((T3,A,T72,A))")& 02078 "HFX_INFO| Interaction Potential: ", "TRUNCATED" 02079 rc_ang = cp_unit_from_cp2k(x_data%potential_parameter%cutoff_radius,"angstrom",error=error) 02080 WRITE(iw,'(T3,A,T61,F20.10)') & 02081 "HFX_INFO| Cutoff Radius [angstrom]: ", rc_ang 02082 CASE(do_hfx_potential_mix_cl_trunc) 02083 WRITE (UNIT=iw,FMT="((T3,A,T65,A))")& 02084 "HFX_INFO| Interaction Potential: ", "TRUNCATED MIX_CL" 02085 rc_ang = cp_unit_from_cp2k(x_data%potential_parameter%cutoff_radius,"angstrom",error=error) 02086 WRITE(iw,'(T3,A,T61,F20.10)') & 02087 "HFX_INFO| Cutoff Radius [angstrom]: ", rc_ang 02088 END SELECT 02089 IF( x_data%periodic_parameter%do_periodic) THEN 02090 IF( x_data%periodic_parameter%mode == -1 ) THEN 02091 WRITE (UNIT=iw,FMT="((T3,A,T77,A))")& 02092 "HFX_INFO| NUMBER_OF_SHELLS: ", "AUTO" 02093 ELSE 02094 WRITE (UNIT=iw,FMT="((T3,A,T61,I20))")& 02095 "HFX_INFO| NUMBER_OF_SHELLS: ", x_data%periodic_parameter%mode 02096 END IF 02097 WRITE (UNIT=iw,FMT="((T3,A,T61,I20))")& 02098 "HFX_INFO| Number of periodic shells considered: ", x_data%periodic_parameter%number_of_shells 02099 WRITE (UNIT=iw,FMT="((T3,A,T61,I20),/)")& 02100 "HFX_INFO| Number of periodic cells considered: ",SIZE(x_data%neighbor_cells) 02101 ELSE 02102 WRITE (UNIT=iw,FMT="((T3,A,T77,A))")& 02103 "HFX_INFO| Number of periodic shells considered: ", "NONE" 02104 WRITE (UNIT=iw,FMT="((T3,A,T77,A),/)")& 02105 "HFX_INFO| Number of periodic cells considered: ", "NONE" 02106 END IF 02107 END IF 02108 CALL cp_print_key_finished_output(iw,logger,hfx_section,& 02109 "HF_INFO", error=error) 02110 END SUBROUTINE hfx_print_info 02111 02112 SUBROUTINE dealloc_containers(x_data, eval_type, error) 02113 TYPE(hfx_type), POINTER :: x_data 02114 INTEGER, INTENT(IN) :: eval_type 02115 TYPE(cp_error_type), INTENT(inout) :: error 02116 02117 CHARACTER(LEN=*), PARAMETER :: routineN = 'dealloc_containers', 02118 routineP = moduleN//':'//routineN 02119 02120 INTEGER :: bin, i, stat 02121 LOGICAL :: failure 02122 02123 failure = .FALSE. 02124 02125 02126 SELECT CASE(eval_type) 02127 CASE (hfx_do_eval_energy) 02128 DO bin=1,SIZE(x_data%maxval_container) 02129 CALL hfx_init_container(x_data%maxval_container(bin), x_data%memory_parameter%actual_memory_usage, & 02130 .FALSE., error) 02131 DEALLOCATE(x_data%maxval_container(bin)%first,STAT=stat) 02132 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02133 END DO 02134 DEALLOCATE(x_data%maxval_container,STAT=stat) 02135 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02136 DEALLOCATE(x_data%maxval_cache,STAT=stat) 02137 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02138 02139 DO bin = 1,SIZE(x_data%integral_containers,2) 02140 DO i=1,64 02141 CALL hfx_init_container(x_data%integral_containers(i,bin), x_data%memory_parameter%actual_memory_usage, & 02142 .FALSE., error) 02143 DEALLOCATE(x_data%integral_containers(i,bin)%first,STAT=stat) 02144 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02145 END DO 02146 END DO 02147 DEALLOCATE(x_data%integral_containers,STAT=stat) 02148 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02149 02150 DEALLOCATE(x_data%integral_caches, STAT=stat) 02151 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02152 CASE (hfx_do_eval_forces) 02153 DO bin=1,SIZE(x_data%maxval_container_forces) 02154 CALL hfx_init_container(x_data%maxval_container_forces(bin), x_data%memory_parameter%actual_memory_usage, & 02155 .FALSE., error) 02156 DEALLOCATE(x_data%maxval_container_forces(bin)%first,STAT=stat) 02157 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02158 END DO 02159 DEALLOCATE(x_data%maxval_container_forces,STAT=stat) 02160 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02161 DEALLOCATE(x_data%maxval_cache_forces,STAT=stat) 02162 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02163 02164 DO bin = 1,SIZE(x_data%integral_containers_forces,2) 02165 DO i=1,64 02166 CALL hfx_init_container(x_data%integral_containers_forces(i,bin), x_data%memory_parameter%actual_memory_usage, & 02167 .FALSE., error) 02168 DEALLOCATE(x_data%integral_containers_forces(i,bin)%first,STAT=stat) 02169 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02170 END DO 02171 END DO 02172 DEALLOCATE(x_data%integral_containers_forces,STAT=stat) 02173 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02174 02175 DEALLOCATE(x_data%integral_caches_forces, STAT=stat) 02176 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02177 02178 END SELECT 02179 02180 END SUBROUTINE dealloc_containers 02181 02182 SUBROUTINE alloc_containers(x_data, bin_size, eval_type, error) 02183 TYPE(hfx_type), POINTER :: x_data 02184 INTEGER, INTENT(IN) :: bin_size, eval_type 02185 TYPE(cp_error_type), INTENT(inout) :: error 02186 02187 CHARACTER(LEN=*), PARAMETER :: routineN = 'alloc_containers', 02188 routineP = moduleN//':'//routineN 02189 02190 INTEGER :: bin, i, stat 02191 LOGICAL :: failure 02192 02193 failure = .FALSE. 02194 02195 02196 SELECT CASE(eval_type) 02197 CASE (hfx_do_eval_energy) 02198 ALLOCATE(x_data%maxval_cache(bin_size),STAT=stat) 02199 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02200 DO bin = 1,bin_size 02201 x_data%maxval_cache(bin)%element_counter = 1 02202 END DO 02203 ALLOCATE(x_data%maxval_container(bin_size),STAT=stat) 02204 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02205 DO bin = 1,bin_size 02206 ALLOCATE(x_data%maxval_container(bin)%first,STAT=stat) 02207 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02208 x_data%maxval_container(bin)%first%prev => NULL() 02209 x_data%maxval_container(bin)%first%next => NULL() 02210 x_data%maxval_container(bin)%current => x_data%maxval_container(bin)%first 02211 x_data%maxval_container(bin)%current%data = 0 02212 x_data%maxval_container(bin)%element_counter = 1 02213 END DO 02214 02215 ALLOCATE(x_data%integral_containers(64,bin_size), STAT=stat) 02216 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02217 ALLOCATE(x_data%integral_caches(64,bin_size), STAT=stat) 02218 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02219 02220 DO bin = 1,bin_size 02221 DO i=1,64 02222 x_data%integral_caches(i,bin)%element_counter = 1 02223 x_data%integral_caches(i,bin)%data = 0 02224 ALLOCATE(x_data%integral_containers(i,bin)%first,STAT=stat) 02225 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02226 x_data%integral_containers(i,bin)%first%prev => NULL() 02227 x_data%integral_containers(i,bin)%first%next => NULL() 02228 x_data%integral_containers(i,bin)%current => x_data%integral_containers(i,bin)%first 02229 x_data%integral_containers(i,bin)%current%data = 0 02230 x_data%integral_containers(i,bin)%element_counter = 1 02231 END DO 02232 END DO 02233 CASE (hfx_do_eval_forces) 02234 ALLOCATE(x_data%maxval_cache_forces(bin_size),STAT=stat) 02235 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02236 DO bin = 1,bin_size 02237 x_data%maxval_cache_forces(bin)%element_counter = 1 02238 END DO 02239 ALLOCATE(x_data%maxval_container_forces(bin_size),STAT=stat) 02240 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02241 DO bin = 1,bin_size 02242 ALLOCATE(x_data%maxval_container_forces(bin)%first,STAT=stat) 02243 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02244 x_data%maxval_container_forces(bin)%first%prev => NULL() 02245 x_data%maxval_container_forces(bin)%first%next => NULL() 02246 x_data%maxval_container_forces(bin)%current => x_data%maxval_container_forces(bin)%first 02247 x_data%maxval_container_forces(bin)%current%data = 0 02248 x_data%maxval_container_forces(bin)%element_counter = 1 02249 END DO 02250 02251 ALLOCATE(x_data%integral_containers_forces(64,bin_size), STAT=stat) 02252 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02253 ALLOCATE(x_data%integral_caches_forces(64,bin_size), STAT=stat) 02254 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02255 02256 DO bin = 1,bin_size 02257 DO i=1,64 02258 x_data%integral_caches_forces(i,bin)%element_counter = 1 02259 x_data%integral_caches_forces(i,bin)%data = 0 02260 ALLOCATE(x_data%integral_containers_forces(i,bin)%first,STAT=stat) 02261 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02262 x_data%integral_containers_forces(i,bin)%first%prev => NULL() 02263 x_data%integral_containers_forces(i,bin)%first%next => NULL() 02264 x_data%integral_containers_forces(i,bin)%current => x_data%integral_containers_forces(i,bin)%first 02265 x_data%integral_containers_forces(i,bin)%current%data = 0 02266 x_data%integral_containers_forces(i,bin)%element_counter = 1 02267 END DO 02268 END DO 02269 02270 END SELECT 02271 02272 END SUBROUTINE alloc_containers 02273 02274 02275 END MODULE hfx_types 02276
1.7.3