|
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_energy_potential 00013 00014 USE atomic_kind_types, ONLY: atomic_kind_type,& 00015 get_atomic_kind_set 00016 USE bibliography, ONLY: cite_reference,& 00017 guidon2008,& 00018 guidon2009 00019 USE cell_types, ONLY: cell_type,& 00020 pbc 00021 USE cp_dbcsr_interface, ONLY: cp_dbcsr_trace 00022 USE cp_dbcsr_types, ONLY: cp_dbcsr_p_type 00023 USE cp_files, ONLY: get_unit_number 00024 USE cp_output_handling, ONLY: cp_p_file,& 00025 cp_print_key_finished_output,& 00026 cp_print_key_should_output,& 00027 cp_print_key_unit_nr 00028 USE cp_para_types, ONLY: cp_para_env_type 00029 USE gamma, ONLY: init_md_ftable 00030 USE hfx_communication, ONLY: distribute_ks_matrix,& 00031 get_atomic_block_maps,& 00032 get_full_density 00033 USE hfx_compression_methods, ONLY: hfx_add_mult_cache_elements,& 00034 hfx_add_single_cache_element,& 00035 hfx_decompress_first_cache,& 00036 hfx_flush_last_cache,& 00037 hfx_get_mult_cache_elements,& 00038 hfx_get_single_cache_element,& 00039 hfx_reset_cache_and_container 00040 USE hfx_contract_block, ONLY: contract_block 00041 USE hfx_libint_interface, ONLY: evaluate_eri 00042 USE hfx_libint_wrapper_types, ONLY: has_iso_c_binding,& 00043 lib_int 00044 USE hfx_load_balance_methods, ONLY: collect_load_balance_info,& 00045 hfx_load_balance,& 00046 hfx_update_load_balance 00047 USE hfx_pair_list_methods, ONLY: build_atomic_pair_list,& 00048 build_pair_list,& 00049 build_pair_list_pgf,& 00050 build_pgf_product_list 00051 USE hfx_screening_methods, ONLY: calc_pair_dist_radii,& 00052 calc_screening_functions,& 00053 update_pmax_mat 00054 USE hfx_types, ONLY: & 00055 alloc_containers, dealloc_containers, hfx_basis_info_type, & 00056 hfx_basis_type, hfx_cache_type, hfx_cell_type, hfx_container_type, & 00057 hfx_create_neighbor_cells, hfx_distribution, hfx_general_type, & 00058 hfx_init_container, hfx_load_balance_type, hfx_memory_type, & 00059 hfx_p_kind, hfx_pgf_list, hfx_pgf_product_list, hfx_potential_type, & 00060 hfx_reset_memory_usage_counter, hfx_screen_coeff_type, & 00061 hfx_screening_type, hfx_task_list_type, hfx_type, init_t_c_g0_lmax, & 00062 log_zero, pair_list_type, pair_set_list_type 00063 USE input_constants, ONLY: do_hfx_potential_mix_cl_trunc,& 00064 do_hfx_potential_truncated,& 00065 hfx_do_eval_energy 00066 USE input_section_types, ONLY: section_vals_type 00067 USE kinds, ONLY: dp,& 00068 int_8 00069 USE machine, ONLY: m_memory,& 00070 m_walltime 00071 USE mathconstants 00072 USE message_passing, ONLY: mp_max,& 00073 mp_sum,& 00074 mp_sync 00075 USE orbital_pointers 00076 USE particle_types, ONLY: particle_type 00077 USE qs_energy_types, ONLY: qs_energy_type 00078 USE qs_environment_types, ONLY: get_qs_env,& 00079 qs_environment_type 00080 USE qs_rho_types, ONLY: qs_rho_type 00081 USE t_c_g0, ONLY: init 00082 USE timings, ONLY: timeset,& 00083 timestop 00084 USE util, ONLY: sort 00085 #include "cp_common_uses.h" 00086 00087 IMPLICIT NONE 00088 PRIVATE 00089 00090 PUBLIC integrate_four_center, coulomb4 00091 00092 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hfx_energy_potential' 00093 00094 !*** 00095 00096 CONTAINS 00097 00098 ! ***************************************************************************** 00120 SUBROUTINE integrate_four_center(qs_env,ks_matrix,energy,rho,hfx_section,para_env,& 00121 geometry_did_change,irep,distribute_fock_matrix,& 00122 ispin, error, do_exx) 00123 00124 TYPE(qs_environment_type), POINTER :: qs_env 00125 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00126 POINTER :: ks_matrix 00127 TYPE(qs_energy_type), POINTER :: energy 00128 TYPE(qs_rho_type), POINTER :: rho 00129 TYPE(section_vals_type), POINTER :: hfx_section 00130 TYPE(cp_para_env_type), POINTER :: para_env 00131 LOGICAL :: geometry_did_change 00132 INTEGER :: irep 00133 LOGICAL, INTENT(IN) :: distribute_fock_matrix 00134 INTEGER, INTENT(IN) :: ispin 00135 TYPE(cp_error_type), INTENT(inout) :: error 00136 LOGICAL, OPTIONAL :: do_exx 00137 00138 CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_four_center', 00139 routineP = moduleN//':'//routineN 00140 00141 INTEGER :: act_atomic_block_offset, act_set_offset, atomic_offset_ac, 00142 atomic_offset_ad, atomic_offset_bc, atomic_offset_bd, bin, 00143 bits_max_val, buffer_left, buffer_size, buffer_start, cache_size, 00144 current_counter, handle, handle_bin, handle_dist_ks, handle_getP, 00145 handle_load, handle_main, i, i_list_ij, i_list_kl, i_set_list_ij, 00146 i_set_list_ij_start, i_set_list_ij_stop, i_set_list_kl, 00147 i_set_list_kl_start, i_set_list_kl_stop, i_thread, iatom, iatom_block, 00148 ikind, iset, iw, j, jatom, jatom_block, jkind, jset, katom, 00149 katom_block, kind_kind_idx, kkind, kset, l_max, latom, latom_block, 00150 lkind, lset, ma, max_am, max_pgf 00151 INTEGER :: max_set, mb, my_bin_size, n_threads, natom, nbits, ncob, 00152 ncos_max, nints, nkind, nneighbors, nseta, nsetb, nsgf_max, nspins, pa, 00153 sgfb, stat, swap_id, unit_id 00154 INTEGER(int_8) :: estimate_to_store_int, mem_compression_counter, 00155 mem_compression_counter_disk, mem_eris, mem_eris_disk, mem_max_val, 00156 my_current_counter, my_istart, n_processes, ncpu, neris_disk, 00157 neris_incore, neris_onthefly, neris_total, shm_mem_compression_counter, 00158 stor_count_int_disk, stor_count_max_val, storage_counter_integrals 00159 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, last_sgf_global 00160 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, lc_max, 00161 lc_min, ld_max, ld_min, npgfa, npgfb, npgfc, npgfd, nsgfa, nsgfb, 00162 nsgfc, nsgfd 00163 INTEGER, DIMENSION(:, :), POINTER :: first_sgfb, nsgfl_a, nsgfl_b, 00164 nsgfl_c, nsgfl_d 00165 LOGICAL :: buffer_overflow, 00166 do_p_screening, do_periodic, 00167 failure, my_do_exx 00168 REAL(dp) :: bintime_start, bintime_stop, cartesian_estimate, 00169 compression_factor, compression_factor_disk, ene_x_aa, ene_x_bb, 00170 eps_schwarz, eps_storage, fac, hf_fraction, max_contraction_val, 00171 max_val1, max_val2, pmax_entry, ra(3), rb(3), rc(3), rd(3), 00172 spherical_estimate, symm_fac 00173 REAL(dp), ALLOCATABLE, DIMENSION(:) :: primitive_integrals 00174 REAL(dp), DIMENSION(:, :), POINTER :: max_contraction, sphi_b, 00175 zeta, zetb, zetc, zetd 00176 TYPE(atomic_kind_type), DIMENSION(:), 00177 POINTER :: atomic_kind_set 00178 TYPE(cp_logger_type), POINTER :: logger 00179 TYPE(hfx_basis_info_type), POINTER :: basis_info 00180 TYPE(hfx_basis_type), DIMENSION(:), 00181 POINTER :: basis_parameter 00182 TYPE(hfx_cache_type), DIMENSION(:), 00183 POINTER :: integral_caches, 00184 integral_caches_disk 00185 TYPE(hfx_cache_type), POINTER :: maxval_cache, 00186 maxval_cache_disk 00187 TYPE(hfx_container_type), DIMENSION(:), 00188 POINTER :: integral_containers, 00189 integral_containers_disk 00190 TYPE(hfx_container_type), POINTER :: maxval_container, 00191 maxval_container_disk 00192 TYPE(hfx_general_type) :: general_parameter 00193 TYPE(hfx_load_balance_type), POINTER :: load_balance_parameter 00194 TYPE(hfx_memory_type), POINTER :: memory_parameter 00195 TYPE(hfx_potential_type) :: potential_parameter 00196 TYPE(hfx_screening_type) :: screening_parameter 00197 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data 00198 TYPE(hfx_type), POINTER :: actual_x_data, 00199 shm_master_x_data 00200 TYPE(particle_type), DIMENSION(:), 00201 POINTER :: particle_set 00202 00203 !$ INTEGER :: omp_get_max_threads,omp_get_thread_num 00204 INTEGER, SAVE :: shm_number_of_p_entries 00205 INTEGER, DIMENSION(:,:), POINTER, SAVE :: shm_is_assoc_atomic_block 00206 INTEGER(int_8) :: shm_neris_total, shm_neris_onthefly, 00207 shm_storage_counter_integrals, shm_neris_incore, 00208 shm_neris_disk, shm_stor_count_int_disk, shm_stor_count_max_val 00209 LOGICAL :: do_disk_storage, use_disk_storage 00210 INTEGER(int_8) :: counter, neris_tmp, nprim_ints, shm_nprim_ints, 00211 mb_size_p, mb_size_f, mb_size_buffers 00212 TYPE(cell_type), POINTER :: cell 00213 CHARACTER(LEN=512) :: error_msg 00214 TYPE(hfx_screen_coeff_type), 00215 DIMENSION(:,:,:,:,:,:), POINTER :: screen_coeffs_pgf, radii_pgf 00216 TYPE(hfx_screen_coeff_type), 00217 DIMENSION(:,:,:,:), POINTER :: screen_coeffs_set 00218 TYPE(hfx_screen_coeff_type), 00219 DIMENSION(:,:), POINTER :: screen_coeffs_kind 00220 REAL(dp) :: rab2, rcd2, 00221 log10_eps_schwarz, 00222 log10_pmax, 00223 ln_10 00224 INTEGER(int_8) :: subtr_size_mb 00225 00226 TYPE(hfx_screen_coeff_type), 00227 DIMENSION(:,:), POINTER :: tmp_R_1, tmp_R_2, tmp_screen_pgf1, tmp_screen_pgf2 00228 00229 INTEGER(int_8) :: memsize_before,memsize_after, 00230 max_val_memory, tmp_i8(8) 00231 REAL(dp) :: max_val2_set 00232 TYPE(lib_int) :: private_lib 00233 REAL(dp), DIMENSION(:), POINTER :: full_density, full_density_beta, full_ks, 00234 full_ks_beta 00235 REAL(dp), DIMENSION(:), ALLOCATABLE :: pbd_buf, pbc_buf, pad_buf, pac_buf, 00236 kbd_buf, kbc_buf, kad_buf, kac_buf 00237 LOGICAL :: treat_lsd_in_core, ks_fully_occ 00238 INTEGER(int_8) :: atom_block, tmp_block, nblocks 00239 TYPE(pair_list_type) :: list_ij, list_kl 00240 REAL(KIND=dp) :: coeffs_kind_max0 00241 TYPE(pair_set_list_type), DIMENSION(:), ALLOCATABLE :: set_list_ij, set_list_kl 00242 TYPE(hfx_p_kind), DIMENSION(:), POINTER :: shm_initial_p 00243 REAL(dp), DIMENSION(:,:), POINTER :: ptr_p_1, ptr_p_2, ptr_p_3, ptr_p_4 00244 INTEGER, DIMENSION(:), POINTER :: shm_block_offset 00245 INTEGER, DIMENSION(:,:), POINTER :: shm_atomic_block_offset 00246 INTEGER, DIMENSION(:,:,:,:), POINTER :: shm_set_offset 00247 INTEGER, DIMENSION(:,:), POINTER :: offset_bd_set, offset_bc_set, offset_ad_set, offset_ac_set 00248 REAL(dp) :: ene_x_aa_diag, ene_x_bb_diag 00249 TYPE(hfx_pgf_list), DIMENSION(:), 00250 ALLOCATABLE :: pgf_list_ij, pgf_list_kl 00251 TYPE(hfx_pgf_product_list), DIMENSION(:), 00252 ALLOCATABLE :: pgf_product_list 00253 REAL(dp), DIMENSION(:,:,:,:), POINTER 00254 :: sphi_a_ext, sphi_b_ext, sphi_c_ext, sphi_d_ext 00255 REAL(dp), DIMENSION(:,:,:), POINTER 00256 :: sphi_a_ext_set, sphi_b_ext_set, sphi_c_ext_set, sphi_d_ext_set 00257 LOGICAL :: do_dynamic_load_balancing, bins_left, do_it 00258 TYPE(hfx_distribution), POINTER :: distribution_energy 00259 INTEGER :: shm_total_bins, shm_task_counter, my_thread_id, my_bin_id , 00260 storage_id 00261 TYPE(hfx_task_list_type), DIMENSION(:), 00262 POINTER :: shm_task_list, tmp_task_list 00263 INTEGER(int_8), DIMENSION(:), 00264 ALLOCATABLE :: tmp_task_list_cost 00265 INTEGER, DIMENSION(:), ALLOCATABLE :: tmp_index 00266 INTEGER(int_8) :: current_ram_counter 00267 REAL(dp), DIMENSION(:), ALLOCATABLE :: ee_work, ee_work2, 00268 ee_buffer1, ee_buffer2, 00269 ee_primitives_tmp 00270 INTEGER, DIMENSION(:), ALLOCATABLE :: nimages 00271 INTEGER :: sphi_d_u1, sphi_d_u2, sphi_d_u3 00272 INTEGER :: sphi_c_u1, sphi_c_u2, sphi_c_u3 00273 INTEGER :: sphi_b_u1, sphi_b_u2, sphi_b_u3 00274 INTEGER :: sphi_a_u1, sphi_a_u2, sphi_a_u3 00275 INTEGER :: iatom_start, iatom_end, jatom_start, jatom_end, katom_start, 00276 katom_end, latom_start, latom_end, my_process_id 00277 REAL(dp), DIMENSION(:,:), POINTER :: shm_pmax_atom, shm_pmax_block 00278 REAL(dp) :: pmax_atom, screen_kind_ij, screen_kind_kl, pmax_blocks 00279 TYPE(cp_dbcsr_p_type),POINTER,DIMENSION(:)::rho_ao 00280 INTEGER::ii,pbc_shells_dummy 00281 REAL(dp), DIMENSION(:), POINTER :: p_work 00282 LOGICAL, DIMENSION(:,:), POINTER :: shm_atomic_pair_list 00283 TYPE(section_vals_type), POINTER :: print_key 00284 LOGICAL :: do_print_load_balance_info 00285 00286 #if defined(__USE_PAT) 00287 INTEGER :: istat 00288 #include "pat_apif.h" 00289 #endif 00290 00291 CALL timeset(routineN,handle) 00292 00293 #if defined(__USE_PAT) 00294 CALL PAT_sampling_state(PAT_STATE_ON,istat) 00295 #endif 00296 00297 rho_ao=>rho%rho_ao 00298 00299 CALL cite_reference(Guidon2008) 00300 CALL cite_reference(Guidon2009) 00301 00302 ! ** This is not very clean, but effective. ispin can only be 2 if we do the beta spin part in core 00303 IF( ispin == 2 ) geometry_did_change = .FALSE. 00304 00305 logger => cp_error_get_logger(error) 00306 00307 IF( geometry_did_change ) THEN 00308 memsize_before=m_memory() 00309 CALL mp_max(memsize_before, para_env%group) 00310 iw = cp_print_key_unit_nr(logger,hfx_section,"HF_INFO",& 00311 extension=".scfLog",error=error) 00312 IF (iw>0) THEN 00313 WRITE (UNIT=iw,FMT="(/,(T3,A,T60,I21))")& 00314 "HFX_MEM_INFO| Est. max. program size before HFX [MB's]:", memsize_before/(1024*1024) 00315 END IF 00316 CALL cp_print_key_finished_output(iw,logger,hfx_section,& 00317 "HF_INFO", error=error) 00318 END IF 00319 00320 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, cell=cell, error=error) 00321 00322 my_do_exx=.FALSE. 00323 IF(PRESENT(do_exx)) my_do_exx=do_exx 00324 NULLIFY(x_data) 00325 IF(my_do_exx) THEN 00326 x_data=>qs_env%mp2_env%ri_rpa%x_data 00327 ELSE 00328 x_data=>qs_env%x_data 00329 END IF 00330 00331 !! Calculate l_max used in fgamma , because init_md_ftable is definitely not thread safe 00332 nkind = SIZE(atomic_kind_set,1) 00333 l_max = 0 00334 DO ikind=1,nkind 00335 l_max = MAX(l_max,MAXVAL(x_data(1,1)%basis_parameter(ikind)%lmax)) 00336 ENDDO 00337 l_max = 4*l_max 00338 CALL init_md_ftable(l_max) 00339 00340 IF(x_data(1,1)%potential_parameter%potential_type == do_hfx_potential_truncated .OR. & 00341 x_data(1,1)%potential_parameter%potential_type == do_hfx_potential_mix_cl_trunc) THEN 00342 IF(l_max>init_t_c_g0_lmax) THEN 00343 IF( para_env%mepos == 0 ) THEN 00344 unit_id = get_unit_number() 00345 OPEN(unit_id,FILE=x_data(1,1)%potential_parameter%filename) 00346 END IF 00347 CALL init(l_max, unit_id, para_env%mepos, para_env%group) 00348 IF(para_env%mepos == 0 ) THEN 00349 CLOSE(unit_id) 00350 END IF 00351 init_t_c_g0_lmax = l_max 00352 END IF 00353 END IF 00354 00355 00356 n_threads = 1 00357 !$ n_threads = omp_get_max_threads() 00358 00359 shm_neris_total = 0 00360 shm_nprim_ints = 0 00361 shm_neris_onthefly = 0 00362 shm_storage_counter_integrals = 0 00363 shm_stor_count_int_disk = 0 00364 shm_neris_incore = 0 00365 shm_neris_disk = 0 00366 shm_stor_count_max_val = 0 00367 00368 !$OMP PARALLEL DEFAULT(PRIVATE) SHARED(qs_env,& 00369 !$OMP x_data,& 00370 !$OMP ks_matrix,& 00371 !$OMP energy,& 00372 !$OMP rho_ao,& 00373 !$OMP hfx_section,& 00374 !$OMP para_env,& 00375 !$OMP geometry_did_change,& 00376 !$OMP irep,& 00377 !$OMP distribute_fock_matrix,& 00378 !$OMP error,& 00379 !$OMP ncoset,& 00380 !$OMP nso,& 00381 !$OMP nco,& 00382 !$OMP full_ks,& 00383 !$OMP full_ks_beta,& 00384 !$OMP n_threads,& 00385 !$OMP full_density,& 00386 !$OMP full_density_beta,& 00387 !$OMP shm_initial_p,& 00388 !$OMP shm_is_assoc_atomic_block,& 00389 !$OMP shm_number_of_p_entries,& 00390 !$OMP shm_neris_total,& 00391 !$OMP shm_neris_onthefly,& 00392 !$OMP shm_storage_counter_integrals,& 00393 !$OMP shm_stor_count_int_disk,& 00394 !$OMP shm_neris_incore,& 00395 !$OMP shm_neris_disk,& 00396 !$OMP shm_nprim_ints,& 00397 !$OMP shm_stor_count_max_val,& 00398 !$OMP cell,& 00399 !$OMP screen_coeffs_set,& 00400 !$OMP screen_coeffs_kind,& 00401 !$OMP screen_coeffs_pgf,& 00402 !$OMP radii_pgf,& 00403 !$OMP nkind,& 00404 !$OMP ispin,& 00405 !$OMP shm_atomic_block_offset,& 00406 !$OMP shm_set_offset,& 00407 !$OMP shm_block_offset,& 00408 !$OMP shm_task_counter,& 00409 !$OMP shm_task_list,& 00410 !$OMP shm_total_bins,& 00411 !$OMP shm_master_x_data,& 00412 !$OMP shm_pmax_atom,& 00413 !$OMP shm_pmax_block,& 00414 !$OMP shm_atomic_pair_list,& 00415 !$OMP shm_mem_compression_counter,& 00416 !$OMP do_print_load_balance_info) 00417 00418 00419 ln_10 = LOG(10.0_dp) 00420 i_thread = 0 00421 !$ i_thread = omp_get_thread_num() 00422 00423 actual_x_data => x_data(irep, i_thread + 1) 00424 !$OMP MASTER 00425 shm_master_x_data => x_data(irep,1) 00426 !$OMP END MASTER 00427 !$OMP BARRIER 00428 00429 00430 00431 do_periodic = actual_x_data%periodic_parameter%do_periodic 00432 00433 IF( do_periodic ) THEN 00434 ! ** Rebuild neighbor lists in case the cell has changed (i.e. NPT MD) 00435 actual_x_data%periodic_parameter%number_of_shells = actual_x_data%periodic_parameter%mode 00436 CALL hfx_create_neighbor_cells(actual_x_data, actual_x_data%periodic_parameter%number_of_shells_from_input,& 00437 cell, i_thread, & 00438 error) 00439 END IF 00440 00441 00442 screening_parameter = actual_x_data%screening_parameter 00443 potential_parameter = actual_x_data%potential_parameter 00444 00445 general_parameter = actual_x_data%general_parameter 00446 load_balance_parameter => actual_x_data%load_balance_parameter 00447 memory_parameter => actual_x_data%memory_parameter 00448 00449 00450 00451 cache_size = memory_parameter%cache_size 00452 bits_max_val = memory_parameter%bits_max_val 00453 00454 basis_parameter => actual_x_data%basis_parameter 00455 basis_info => actual_x_data%basis_info 00456 00457 treat_lsd_in_core = general_parameter%treat_lsd_in_core 00458 00459 ncpu = para_env%num_pe 00460 n_processes = ncpu * n_threads 00461 00462 !! initalize some counters 00463 neris_total = 0_int_8 00464 neris_incore = 0_int_8 00465 neris_disk = 0_int_8 00466 neris_onthefly = 0_int_8 00467 mem_eris = 0_int_8 00468 mem_eris_disk = 0_int_8 00469 mem_max_val = 0_int_8 00470 compression_factor = 0.0_dp 00471 compression_factor_disk = 0.0_dp 00472 nprim_ints = 0_int_8 00473 neris_tmp = 0_int_8 00474 max_val_memory = 1_int_8 00475 00476 max_am = basis_info%max_am 00477 00478 do_p_screening = screening_parameter%do_initial_p_screening 00479 max_set = basis_info%max_set 00480 CALL get_qs_env(qs_env=qs_env,& 00481 atomic_kind_set=atomic_kind_set,& 00482 particle_set=particle_set,& 00483 error=error) 00484 00485 natom = SIZE(particle_set,1) 00486 00487 ALLOCATE(kind_of(natom),STAT=stat) 00488 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00489 00490 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 00491 kind_of=kind_of) 00492 00493 !! precompute maximum nco and allocate scratch 00494 ncos_max=0 00495 nsgf_max=0 00496 DO iatom=1,natom 00497 ikind = kind_of(iatom) 00498 nseta = basis_parameter(ikind)%nset 00499 npgfa => basis_parameter(ikind)%npgf 00500 la_max => basis_parameter(ikind)%lmax 00501 nsgfa => basis_parameter(ikind)%nsgf 00502 DO iset = 1, nseta 00503 ncos_max = MAX(ncos_max,ncoset(la_max(iset))) 00504 nsgf_max = MAX(nsgf_max,nsgfa(iset)) 00505 ENDDO 00506 ENDDO 00507 !! Allocate the arrays for the integrals. 00508 ALLOCATE(primitive_integrals(nsgf_max**4),STAT=stat) 00509 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00510 primitive_integrals = 0.0_dp 00511 00512 ALLOCATE(pbd_buf(nsgf_max**2),STAT=stat) 00513 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00514 ALLOCATE(pbc_buf(nsgf_max**2),STAT=stat) 00515 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00516 ALLOCATE(pad_buf(nsgf_max**2),STAT=stat) 00517 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00518 ALLOCATE(pac_buf(nsgf_max**2),STAT=stat) 00519 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00520 ALLOCATE(kbd_buf(nsgf_max**2),STAT=stat) 00521 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00522 ALLOCATE(kbc_buf(nsgf_max**2),STAT=stat) 00523 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00524 ALLOCATE(kad_buf(nsgf_max**2),STAT=stat) 00525 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00526 ALLOCATE(kac_buf(nsgf_max**2),STAT=stat) 00527 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00528 ALLOCATE(ee_work(ncos_max**4),STAT=stat) 00529 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00530 ALLOCATE(ee_work2(ncos_max**4),STAT=stat) 00531 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00532 ALLOCATE(ee_buffer1(ncos_max**4),STAT=stat) 00533 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00534 ALLOCATE(ee_buffer2(ncos_max**4),STAT=stat) 00535 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00536 ALLOCATE(ee_primitives_tmp(nsgf_max**4),STAT=stat) ! XXXXX could be wrong 00537 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00538 00539 IF( .NOT. has_iso_c_binding) THEN 00540 ALLOCATE(p_work(nco(max_am)**4), STAT=stat) 00541 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00542 END IF 00543 00544 nspins = qs_env%dft_control%nspins 00545 00546 ALLOCATE(max_contraction(max_set,natom),STAT=stat) 00547 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00548 00549 max_contraction=0.0_dp 00550 max_pgf = 0 00551 DO jatom=1,natom 00552 jkind = kind_of(jatom) 00553 lb_max => basis_parameter(jkind)%lmax 00554 nsetb = basis_parameter(jkind)%nset 00555 npgfb => basis_parameter(jkind)%npgf 00556 first_sgfb => basis_parameter(jkind)%first_sgf 00557 sphi_b => basis_parameter(jkind)%sphi 00558 nsgfb => basis_parameter(jkind)%nsgf 00559 DO jset = 1,nsetb 00560 ! takes the primitive to contracted transformation into account 00561 ncob = npgfb(jset)*ncoset(lb_max(jset)) 00562 sgfb = first_sgfb(1,jset) 00563 ! if the primitives are assumed to be all of max_val2, max_val2*p2s_b becomes 00564 ! the maximum value after multiplication with sphi_b 00565 max_contraction(jset,jatom) = MAXVAL((/(SUM(ABS(sphi_b(1:ncob,i))),i=sgfb,sgfb+nsgfb(jset)-1)/)) 00566 max_pgf = MAX(max_pgf,npgfb(jset)) 00567 ENDDO 00568 ENDDO 00569 00570 ! ** Allocate buffers for pgf_lists 00571 nneighbors = SIZE(actual_x_data%neighbor_cells) 00572 ALLOCATE(pgf_list_ij(max_pgf**2), STAT=stat) 00573 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00574 ALLOCATE(pgf_list_kl(max_pgf**2), STAT=stat) 00575 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00576 ALLOCATE(pgf_product_list(nneighbors**3), STAT=stat) 00577 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00578 ALLOCATE(nimages(max_pgf**2),STAT=stat) 00579 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00580 00581 DO i=1,max_pgf**2 00582 ALLOCATE(pgf_list_ij(i)%image_list(nneighbors), STAT=stat) 00583 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00584 ALLOCATE(pgf_list_kl(i)%image_list(nneighbors), STAT=stat) 00585 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00586 END DO 00587 !$OMP BARRIER 00588 !$OMP MASTER 00589 !! Calculate helper array that stores if a certain atomic pair is associated in the KS matrix 00590 IF( geometry_did_change ) THEN 00591 CALL get_atomic_block_maps(ks_matrix(1)%matrix,basis_parameter,kind_of,& 00592 shm_master_x_data%is_assoc_atomic_block,& 00593 shm_master_x_data%number_of_p_entries ,& 00594 para_env, & 00595 shm_master_x_data%atomic_block_offset,& 00596 shm_master_x_data%set_offset,& 00597 shm_master_x_data%block_offset,& 00598 shm_master_x_data%map_atoms_to_cpus,& 00599 nkind,& 00600 error) 00601 00602 shm_is_assoc_atomic_block => shm_master_x_data%is_assoc_atomic_block 00603 00604 !! Get occupation of KS-matrix 00605 ks_fully_occ = .TRUE. 00606 outer: DO iatom=1,natom 00607 DO jatom=iatom,natom 00608 IF(shm_is_assoc_atomic_block(jatom,iatom) == 0 ) THEN 00609 ks_fully_occ = .FALSE. 00610 EXIT outer 00611 END IF 00612 END DO 00613 END DO outer 00614 00615 IF ( .NOT. ks_fully_occ) THEN 00616 WRITE(error_msg,'(A,A,A)') "The Kohn Sham matrix is not 100% occupied. This "//& 00617 "may result in uncorrect Hartree-Fock results. Try to decrease EPS_PGF_ORB "//& 00618 "and EPS_FILTER_MATRIX in the QS section. " 00619 CALL cp_assert(.FALSE.,cp_warning_level,cp_assertion_failed,routineP,& 00620 error_msg,& 00621 error,& 00622 only_ionode=.TRUE.) 00623 END IF 00624 END IF 00625 00626 ! ** Set pointers 00627 shm_number_of_p_entries = shm_master_x_data%number_of_p_entries 00628 shm_is_assoc_atomic_block => shm_master_x_data%is_assoc_atomic_block 00629 shm_atomic_block_offset => shm_master_x_data%atomic_block_offset 00630 shm_set_offset => shm_master_x_data%set_offset 00631 shm_block_offset => shm_master_x_data%block_offset 00632 !$OMP END MASTER 00633 !$OMP BARRIER 00634 00635 ! ** Reset storage counter given by MAX_MEMORY by subtracting all buffers 00636 ! ** Fock and density Matrices (shared) 00637 subtr_size_mb = 2_int_8*shm_block_offset(ncpu+1) 00638 ! ** if non restricted ==> alpha/beta spin 00639 IF (.NOT. treat_lsd_in_core ) THEN 00640 IF( nspins == 2) subtr_size_mb = subtr_size_mb * 2_int_8 00641 END IF 00642 ! ** Initial P only MAX(alpha,beta) (shared) 00643 IF(do_p_screening .OR. screening_parameter%do_p_screening_forces) THEN 00644 subtr_size_mb = subtr_size_mb + memory_parameter%size_p_screen 00645 END IF 00646 ! ** In core forces require their own initial P 00647 IF(screening_parameter%do_p_screening_forces) THEN 00648 IF( memory_parameter%treat_forces_in_core) THEN 00649 subtr_size_mb = subtr_size_mb + memory_parameter%size_p_screen 00650 END IF 00651 END IF 00652 ! ** primitive buffer (not shared by the threads) 00653 subtr_size_mb = subtr_size_mb + nsgf_max**4*n_threads 00654 ! ** density + fock buffers 00655 subtr_size_mb = subtr_size_mb + 8_int_8*nsgf_max**2*n_threads 00656 ! ** screening functions (shared) 00657 ! ** coeffs_pgf 00658 subtr_size_mb = subtr_size_mb + max_pgf**2*max_set**2*nkind**2 00659 ! ** coeffs_set 00660 subtr_size_mb = subtr_size_mb + max_set**2*nkind**2 00661 ! ** coeffs_kind 00662 subtr_size_mb = subtr_size_mb + nkind**2 00663 ! ** radii_pgf 00664 subtr_size_mb = subtr_size_mb + max_pgf**2*max_set**2*nkind**2 00665 00666 ! ** is_assoc (shared) 00667 subtr_size_mb = subtr_size_mb + natom**2 00668 00669 ! ** pmax_atom (shared) 00670 IF( do_p_screening) THEN 00671 subtr_size_mb = subtr_size_mb + natom**2 00672 END IF 00673 IF(screening_parameter%do_p_screening_forces) THEN 00674 IF( memory_parameter%treat_forces_in_core) THEN 00675 subtr_size_mb = subtr_size_mb + natom**2 00676 END IF 00677 END IF 00678 00679 ! ** Convert into MB's 00680 subtr_size_mb = subtr_size_mb * 8_int_8/1024_int_8/1024_int_8 00681 00682 ! ** Subtracting all these buffers from from MAX_MEMORY yields the amount 00683 ! ** of RAM that is left for the compressed integrals. When using threads 00684 ! ** all the available memory is shared among all n_threads. i.e. the faster 00685 ! ** ones can steal from the slower ones 00686 00687 CALL hfx_reset_memory_usage_counter(memory_parameter, subtr_size_mb) 00688 00689 use_disk_storage = .FALSE. 00690 counter = 0_int_8 00691 do_disk_storage = memory_parameter%do_disk_storage 00692 IF (do_disk_storage ) THEN 00693 maxval_container_disk => actual_x_data%maxval_container_disk 00694 maxval_cache_disk => actual_x_data%maxval_cache_disk 00695 00696 integral_containers_disk => actual_x_data%integral_containers_disk 00697 integral_caches_disk => actual_x_data%integral_caches_disk 00698 END IF 00699 00700 00701 IF(geometry_did_change) THEN 00702 memory_parameter%ram_counter = HUGE(memory_parameter%ram_counter) 00703 END IF 00704 00705 IF( geometry_did_change) THEN 00706 memory_parameter%recalc_forces = .TRUE. 00707 ELSE 00708 IF( .NOT. memory_parameter%treat_forces_in_core ) memory_parameter%recalc_forces = .TRUE. 00709 END IF 00710 00711 !! Get screening parameter 00712 eps_schwarz = screening_parameter%eps_schwarz 00713 IF( eps_schwarz <= 0.0_dp) THEN 00714 log10_eps_schwarz = log_zero 00715 ELSE 00716 log10_eps_schwarz = LOG10(eps_schwarz) 00717 END IF 00718 !! get storage epsilon 00719 eps_storage = eps_schwarz*memory_parameter%eps_storage_scaling 00720 00721 00722 !! If we have a hybrid functional, we may need only a fraction of exact exchange 00723 hf_fraction = general_parameter%fraction 00724 00725 !! The number of integrals that fit into the given MAX_MEMORY 00726 00727 !! Parameters related to the potential 1/r, erf(wr)/r, erfc(wr/r) 00728 potential_parameter = actual_x_data%potential_parameter 00729 00730 !! Variable to check if we calculate the integrals in-core or on the fly 00731 !! If TRUE -> on the fly 00732 IF(.NOT. memory_parameter%do_all_on_the_fly) THEN 00733 buffer_overflow = .FALSE. 00734 ELSE 00735 buffer_overflow = .TRUE. 00736 END IF 00737 logger => cp_error_get_logger(error) 00738 00739 00740 private_lib = actual_x_data%lib 00741 00742 !! Helper array to map local basis function indeces to global ones 00743 ALLOCATE(last_sgf_global(0:natom),STAT=stat) 00744 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00745 last_sgf_global(0)=0 00746 DO iatom=1,natom 00747 ikind = kind_of(iatom) 00748 last_sgf_global(iatom) = last_sgf_global(iatom-1)+ basis_parameter(ikind)%nsgf_total 00749 END DO 00750 !$OMP BARRIER 00751 !$OMP MASTER 00752 !! Let master thread get the density (avoid problems with MPI) 00753 !! Get the full density from all the processors 00754 NULLIFY(full_density) 00755 NULLIFY(full_density_beta) 00756 ALLOCATE(full_density(shm_block_offset(ncpu+1)),STAT=stat) 00757 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00758 IF( .NOT. treat_lsd_in_core .OR. nspins == 1 ) THEN 00759 CALL timeset(routineN//"_getP",handle_getP) 00760 CALL get_full_density(para_env, full_density, rho_ao(ispin)%matrix, shm_number_of_p_entries,& 00761 shm_master_x_data%block_offset, natom, & 00762 kind_of, basis_parameter, get_max_vals_spin=.FALSE., error=error) 00763 00764 IF(nspins == 2) THEN 00765 ALLOCATE(full_density_beta(shm_block_offset(ncpu+1)),STAT=stat) 00766 CALL get_full_density(para_env, full_density_beta, rho_ao(2)%matrix, shm_number_of_p_entries,& 00767 shm_master_x_data%block_offset, natom, & 00768 kind_of, basis_parameter, get_max_vals_spin=.FALSE., error=error) 00769 END IF 00770 CALL timestop(handle_getP) 00771 00772 !! Calculate the max values of the density matrix actual_pmax stores the data from the actual density matrix 00773 !! and x_data%initial_p stores the same for the initial guess. The initial guess is updated only in the case of 00774 !! a changed geometry 00775 NULLIFY(shm_initial_p) 00776 IF( do_p_screening) THEN 00777 shm_initial_p => shm_master_x_data%initial_p 00778 shm_pmax_atom => shm_master_x_data%pmax_atom 00779 IF(geometry_did_change ) THEN 00780 CALL update_pmax_mat(shm_master_x_data%initial_p, & 00781 shm_master_x_data%map_atom_to_kind_atom, & 00782 shm_master_x_data%set_offset,& 00783 shm_master_x_data%atomic_block_offset,& 00784 shm_pmax_atom,& 00785 full_density,full_density_beta,& 00786 natom, kind_of, basis_parameter, & 00787 nkind, shm_master_x_data%is_assoc_atomic_block) 00788 END IF 00789 END IF 00790 ELSE 00791 IF( do_p_screening ) THEN 00792 CALL timeset(routineN//"_getP",handle_getP) 00793 CALL get_full_density(para_env, full_density, rho_ao(1)%matrix, shm_number_of_p_entries,& 00794 shm_master_x_data%block_offset, natom, & 00795 kind_of, basis_parameter, get_max_vals_spin=.TRUE., & 00796 rho_beta=rho_ao(2)%matrix, error=error) 00797 CALL timestop(handle_getP) 00798 00799 !! Calculate the max values of the density matrix actual_pmax stores the data from the actual density matrix 00800 !! and x_data%initial_p stores the same for the initial guess. The initial guess is updated only in the case of 00801 !! a changed geometry 00802 NULLIFY(shm_initial_p) 00803 shm_initial_p => actual_x_data%initial_p 00804 shm_pmax_atom => shm_master_x_data%pmax_atom 00805 IF(geometry_did_change ) THEN 00806 CALL update_pmax_mat(shm_master_x_data%initial_p, & 00807 shm_master_x_data%map_atom_to_kind_atom, & 00808 shm_master_x_data%set_offset,& 00809 shm_master_x_data%atomic_block_offset,& 00810 shm_pmax_atom,& 00811 full_density,full_density_beta,& 00812 natom, kind_of, basis_parameter, & 00813 nkind, shm_master_x_data%is_assoc_atomic_block) 00814 END IF 00815 END IF 00816 ! ** Now get the density(ispin) 00817 CALL get_full_density(para_env, full_density, rho_ao(ispin)%matrix, shm_number_of_p_entries,& 00818 shm_master_x_data%block_offset, natom, & 00819 kind_of, basis_parameter, get_max_vals_spin=.FALSE., error=error) 00820 END IF 00821 00822 NULLIFY(full_ks, full_ks_beta) 00823 ALLOCATE(shm_master_x_data%full_ks_alpha(shm_block_offset(ncpu+1)),STAT=stat) 00824 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00825 full_ks => shm_master_x_data%full_ks_alpha 00826 full_ks = 0.0_dp 00827 00828 00829 IF (.NOT. treat_lsd_in_core ) THEN 00830 IF(nspins==2) THEN 00831 ALLOCATE(shm_master_x_data%full_ks_beta(shm_block_offset(ncpu+1)),STAT=stat) 00832 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00833 full_ks_beta => shm_master_x_data%full_ks_beta 00834 full_ks_beta = 0.0_dp 00835 END IF 00836 END IF 00837 00838 !! Initialize schwarz screening matrices for near field estimates and boxing screening matrices 00839 !! for far field estimates. The update is only performed if the geomtry of the system changed. 00840 !! If the system is periodic, then the corresponding routines are called and some variables 00841 !! are initialized 00842 00843 !$OMP END MASTER 00844 !$OMP BARRIER 00845 00846 00847 IF( .NOT. shm_master_x_data%screen_funct_is_initialized ) THEN 00848 CALL calc_pair_dist_radii(qs_env, basis_parameter,& 00849 shm_master_x_data%pair_dist_radii_pgf, max_set, max_pgf, eps_schwarz,& 00850 n_threads, i_thread, error) 00851 !$OMP BARRIER 00852 CALL calc_screening_functions(qs_env, basis_parameter, private_lib, shm_master_x_data%potential_parameter,& 00853 shm_master_x_data%screen_funct_coeffs_set,& 00854 shm_master_x_data%screen_funct_coeffs_kind, & 00855 shm_master_x_data%screen_funct_coeffs_pgf, & 00856 shm_master_x_data%pair_dist_radii_pgf,& 00857 max_set, max_pgf, n_threads, i_thread, p_work, error) 00858 00859 !$OMP MASTER 00860 shm_master_x_data%screen_funct_is_initialized = .TRUE. 00861 !$OMP END MASTER 00862 END IF 00863 !$OMP BARRIER 00864 00865 !$OMP MASTER 00866 screen_coeffs_set => shm_master_x_data%screen_funct_coeffs_set 00867 screen_coeffs_kind => shm_master_x_data%screen_funct_coeffs_kind 00868 screen_coeffs_pgf => shm_master_x_data%screen_funct_coeffs_pgf 00869 radii_pgf => shm_master_x_data%pair_dist_radii_pgf 00870 !$OMP END MASTER 00871 !$OMP BARRIER 00872 00873 !! Initialize a prefactor depending on the fraction and the number of spins 00874 IF(nspins == 1) THEN 00875 fac = 0.5_dp * hf_fraction 00876 ELSE 00877 fac = 1.0_dp * hf_fraction 00878 END IF 00879 00880 !! Call routines that distribute the load on all processes. If we want to screen on a initial density matrix, there is 00881 !! an optional parameter. Of course, this is only done if the geometry did change 00882 !$OMP BARRIER 00883 !$OMP MASTER 00884 CALL timeset(routineN//"_load",handle_load) 00885 !$OMP END MASTER 00886 !$OMP BARRIER 00887 IF( geometry_did_change ) THEN 00888 IF( actual_x_data%b_first_load_balance_energy ) THEN 00889 CALL hfx_load_balance(actual_x_data,eps_schwarz,particle_set,max_set,para_env,& 00890 screen_coeffs_set,screen_coeffs_kind,& 00891 shm_is_assoc_atomic_block,do_periodic,load_balance_parameter, & 00892 kind_of, basis_parameter, shm_initial_p, shm_pmax_atom, i_thread, n_threads, & 00893 cell, do_p_screening, actual_x_data%map_atom_to_kind_atom, & 00894 nkind, hfx_do_eval_energy, shm_pmax_block, use_virial=.FALSE., error=error) 00895 actual_x_data%b_first_load_balance_energy = .FALSE. 00896 ELSE 00897 CALL hfx_update_load_balance(actual_x_data,para_env, & 00898 load_balance_parameter, & 00899 i_thread, n_threads, hfx_do_eval_energy, & 00900 error) 00901 END IF 00902 END IF 00903 !$OMP BARRIER 00904 !$OMP MASTER 00905 CALL timestop(handle_load) 00906 !$OMP END MASTER 00907 !$OMP BARRIER 00908 00909 !! Start caluclating integrals of the form (ab|cd) or (ij|kl) 00910 !! In order to do so, there is a main four-loop structre that takes into account the two symmetries 00911 !! 00912 !! (ab|cd) = (ba|cd) = (ab|dc) = (ba|dc) 00913 !! 00914 !! by iterating in the following way 00915 !! 00916 !! DO iatom=1,natom and DO katom=1,natom 00917 !! DO jatom=iatom,natom DO latom=katom,natom 00918 !! 00919 !! The third symmetry 00920 !! 00921 !! (ab|cd) = (cd|ab) 00922 !! 00923 !! is taken into account by the following criterion: 00924 !! 00925 !! IF(katom+latom<=iatom+jatom) THEN 00926 !! IF( ((iatom+jatom).EQ.(katom+latom) ) .AND.(katom<iatom)) CYCLE 00927 !! 00928 !! Depending on the degeneracy of an integral the exchange contribution is multiplied by a corresponding 00929 !! factor ( symm_fac ). 00930 !! 00931 !! If one quartet does not pass the screening we CYCLE on the outer most possible loop. Thats why we use 00932 !! different hierarchies of short range screening matrices. 00933 !! 00934 !! If we do a parallel run, each process owns a unique array of workloads. Here, a workload is 00935 !! defined as : 00936 !! 00937 !! istart, jstart, kstart, lstart, number_of_atom_quartets, initial_cpu_id 00938 !! 00939 !! This tells the process where to start the main loops and how many bunches of integrals it has to 00940 !! calculate. The original parallelization is a simple modulo distribution that is binned and 00941 !! optimized in the load_balance routines. Since the Monte Carlo routines can swap processors, 00942 !! we need to know which was the inital cpu_id. 00943 !! Furthermore, the indices iatom, jatom, katom, latom have to be set to istart, jstart, kstart and 00944 !! lstart only the first time the loop is executed. All subsequent loops have to start with one or 00945 !! iatom and katom respectively. Therefore, we use flags like first_j_loop etc. 00946 00947 00948 do_dynamic_load_balancing = .TRUE. 00949 00950 IF( n_threads == 1 .OR. do_disk_storage ) do_dynamic_load_balancing = .FALSE. 00951 00952 IF( do_dynamic_load_balancing ) THEN 00953 my_bin_size = SIZE(actual_x_data%distribution_energy) 00954 ELSE 00955 my_bin_size = 1 00956 END IF 00957 !! We do not need the containers if MAX_MEM = 0 00958 IF(.NOT. memory_parameter%do_all_on_the_fly) THEN 00959 !! IF new md step -> reinitialize containers 00960 IF(geometry_did_change) THEN 00961 CALL dealloc_containers(actual_x_data, hfx_do_eval_energy, error) 00962 CALL alloc_containers(actual_x_data, my_bin_size, hfx_do_eval_energy, error) 00963 00964 DO bin = 1,my_bin_size 00965 maxval_container => actual_x_data%maxval_container(bin) 00966 integral_containers => actual_x_data%integral_containers(:,bin) 00967 CALL hfx_init_container(maxval_container, memory_parameter%actual_memory_usage, .FALSE., error) 00968 DO i=1,64 00969 CALL hfx_init_container(integral_containers(i),memory_parameter%actual_memory_usage, .FALSE., error) 00970 END DO 00971 END DO 00972 END IF 00973 00974 !! Decompress the first cache for maxvals and integrals 00975 IF( .NOT. geometry_did_change) THEN 00976 DO bin = 1,my_bin_size 00977 maxval_cache => actual_x_data%maxval_cache(bin) 00978 maxval_container => actual_x_data%maxval_container(bin) 00979 integral_caches => actual_x_data%integral_caches(:,bin) 00980 integral_containers => actual_x_data%integral_containers(:,bin) 00981 CALL hfx_decompress_first_cache(bits_max_val, maxval_cache, maxval_container, & 00982 memory_parameter%actual_memory_usage, .FALSE.) 00983 DO i=1,64 00984 CALL hfx_decompress_first_cache(i, integral_caches(i), integral_containers(i), & 00985 memory_parameter%actual_memory_usage,.FALSE.) 00986 END DO 00987 END DO 00988 END IF 00989 END IF 00990 00991 !! Since the I/O routines are no thread-safe, i.e. the procedure to get the unit number, put a lock here 00992 !$OMP CRITICAL 00993 !! If we do disk storage, we have to initialize the containers/caches 00994 IF(do_disk_storage) THEN 00995 !! IF new md step -> reinitialize containers 00996 IF(geometry_did_change) THEN 00997 CALL hfx_init_container(maxval_container_disk, memory_parameter%actual_memory_usage_disk, do_disk_storage, error) 00998 DO i=1,64 00999 CALL hfx_init_container(integral_containers_disk(i),memory_parameter%actual_memory_usage_disk, do_disk_storage, error) 01000 END DO 01001 END IF 01002 !! Decompress the first cache for maxvals and integrals 01003 IF( .NOT. geometry_did_change) THEN 01004 CALL hfx_decompress_first_cache(bits_max_val, maxval_cache_disk, maxval_container_disk, & 01005 memory_parameter%actual_memory_usage_disk, .TRUE.) 01006 DO i=1,64 01007 CALL hfx_decompress_first_cache(i, integral_caches_disk(i), integral_containers_disk(i), & 01008 memory_parameter%actual_memory_usage_disk, .TRUE.) 01009 END DO 01010 END IF 01011 END IF 01012 !$OMP END CRITICAL 01013 01014 !$OMP BARRIER 01015 !$OMP MASTER 01016 01017 IF( do_dynamic_load_balancing ) THEN 01018 ! ** Lets contstruct the task list 01019 shm_total_bins = 0 01020 DO i = 1,n_threads 01021 shm_total_bins = shm_total_bins + SIZE(x_data(irep, i)%distribution_energy) 01022 END DO 01023 ALLOCATE(tmp_task_list(shm_total_bins), STAT=stat) 01024 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01025 shm_task_counter = 0 01026 DO i = 1,n_threads 01027 DO bin = 1,SIZE(x_data(irep,i)%distribution_energy) 01028 shm_task_counter = shm_task_counter + 1 01029 tmp_task_list(shm_task_counter)%thread_id = i 01030 tmp_task_list(shm_task_counter)%bin_id = bin 01031 tmp_task_list(shm_task_counter)%cost = x_data(irep,i)%distribution_energy(bin)%cost 01032 END DO 01033 END DO 01034 01035 ! ** Now sort the task list 01036 ALLOCATE(tmp_task_list_cost(shm_total_bins),STAT=stat) 01037 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01038 ALLOCATE(tmp_index(shm_total_bins),STAT=stat) 01039 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01040 01041 DO i = 1,shm_total_bins 01042 tmp_task_list_cost(i) = tmp_task_list(i)%cost 01043 END DO 01044 01045 CALL sort(tmp_task_list_cost,shm_total_bins,tmp_index) 01046 01047 ALLOCATE(shm_master_x_data%task_list(shm_total_bins), STAT=stat) 01048 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01049 01050 DO i = 1,shm_total_bins 01051 shm_master_x_data%task_list(i) = tmp_task_list(tmp_index(shm_total_bins - i +1)) 01052 END DO 01053 01054 shm_task_list => shm_master_x_data%task_list 01055 shm_task_counter = 0 01056 01057 DEALLOCATE(tmp_task_list_cost, tmp_index, tmp_task_list, STAT=stat) 01058 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01059 END IF 01060 !$OMP END MASTER 01061 !$OMP BARRIER 01062 01063 IF( my_bin_size > 0 ) THEN 01064 maxval_container => actual_x_data%maxval_container(1) 01065 maxval_cache => actual_x_data%maxval_cache(1) 01066 01067 integral_containers => actual_x_data%integral_containers(:,1) 01068 integral_caches => actual_x_data%integral_caches(:,1) 01069 END IF 01070 01071 !$OMP BARRIER 01072 !$OMP MASTER 01073 CALL timeset(routineN//"_main",handle_main) 01074 !$OMP END MASTER 01075 !$OMP BARRIER 01076 01077 01078 coeffs_kind_max0=MAXVAL(screen_coeffs_kind(:,:)%x(2)) 01079 ALLOCATE(set_list_ij((max_set*load_balance_parameter%block_size)**2),STAT=stat) 01080 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01081 ALLOCATE(set_list_kl((max_set*load_balance_parameter%block_size)**2),STAT=stat) 01082 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01083 01084 !$OMP BARRIER 01085 !$OMP MASTER 01086 01087 !! precalculate maximum density matrix elements in blocks 01088 actual_x_data%pmax_block = 0.0_dp 01089 shm_pmax_block => actual_x_data%pmax_block 01090 IF( do_p_screening ) THEN 01091 DO iatom_block = 1,SIZE(actual_x_data%blocks) 01092 iatom_start = actual_x_data%blocks(iatom_block)%istart 01093 iatom_end = actual_x_data%blocks(iatom_block)%iend 01094 DO jatom_block = 1,SIZE(actual_x_data%blocks) 01095 jatom_start = actual_x_data%blocks(jatom_block)%istart 01096 jatom_end = actual_x_data%blocks(jatom_block)%iend 01097 shm_pmax_block(iatom_block,jatom_block) = MAXVAL(shm_pmax_atom(iatom_start:iatom_end,jatom_start:jatom_end)) 01098 END DO 01099 END DO 01100 END IF 01101 shm_atomic_pair_list => actual_x_data%atomic_pair_list 01102 IF( geometry_did_change) THEN 01103 CALL build_atomic_pair_list(natom, shm_atomic_pair_list, kind_of,basis_parameter,particle_set, & 01104 do_periodic,screen_coeffs_kind,coeffs_kind_max0,log10_eps_schwarz,cell,& 01105 actual_x_data%blocks) 01106 END IF 01107 01108 my_bin_size = SIZE(actual_x_data%distribution_energy) 01109 ! reset timings for the new SCF round 01110 IF (geometry_did_change) THEN 01111 DO bin=1,my_bin_size 01112 actual_x_data%distribution_energy(bin)%time_first_scf=0.0_dp 01113 actual_x_data%distribution_energy(bin)%time_other_scf=0.0_dp 01114 actual_x_data%distribution_energy(bin)%time_forces=0.0_dp 01115 END DO 01116 ENDIF 01117 !$OMP END MASTER 01118 !$OMP BARRIER 01119 01120 my_bin_size = SIZE(actual_x_data%distribution_energy) 01121 nblocks = load_balance_parameter%nblocks 01122 01123 bins_left = .TRUE. 01124 do_it = .TRUE. 01125 bin = 0 01126 DO WHILE (bins_left) 01127 IF( .NOT. do_dynamic_load_balancing) THEN 01128 bin = bin + 1 01129 IF( bin > my_bin_size ) THEN 01130 do_it = .FALSE. 01131 bins_left = .FALSE. 01132 ELSE 01133 do_it = .TRUE. 01134 bins_left = .TRUE. 01135 distribution_energy => actual_x_data%distribution_energy(bin) 01136 END IF 01137 ELSE 01138 !$OMP CRITICAL 01139 shm_task_counter = shm_task_counter + 1 01140 IF( shm_task_counter <= shm_total_bins) THEN 01141 my_thread_id = shm_task_list(shm_task_counter)%thread_id 01142 my_bin_id = shm_task_list(shm_task_counter)%bin_id 01143 maxval_cache => x_data(irep,my_thread_id)%maxval_cache(my_bin_id) 01144 maxval_container => x_data(irep,my_thread_id)%maxval_container(my_bin_id) 01145 integral_caches => x_data(irep,my_thread_id)%integral_caches(:,my_bin_id) 01146 integral_containers => x_data(irep,my_thread_id)%integral_containers(:,my_bin_id) 01147 distribution_energy => x_data(irep,my_thread_id)%distribution_energy(my_bin_id) 01148 do_it = .TRUE. 01149 bins_left = .TRUE. 01150 IF( geometry_did_change) THEN 01151 distribution_energy%ram_counter = HUGE(distribution_energy%ram_counter) 01152 END IF 01153 counter = 0_Int_8 01154 ELSE 01155 do_it = .FALSE. 01156 bins_left = .FALSE. 01157 END IF 01158 !$OMP END CRITICAL 01159 END IF 01160 01161 IF( .NOT. do_it ) CYCLE 01162 !$OMP MASTER 01163 CALL timeset(routineN//"_bin",handle_bin) 01164 !$OMP END MASTER 01165 bintime_start=m_walltime() 01166 my_istart = distribution_energy%istart 01167 my_current_counter = 0 01168 IF(distribution_energy%number_of_atom_quartets == 0 .OR. & 01169 my_istart == -1_int_8) my_istart = nblocks **4 01170 atomic_blocks: DO atom_block = my_istart,nblocks**4-1,n_processes 01171 latom_block = MODULO(atom_block,nblocks)+1 01172 tmp_block = atom_block/nblocks 01173 katom_block = MODULO(tmp_block,nblocks)+1 01174 IF(latom_block<katom_block) CYCLE 01175 tmp_block = tmp_block/nblocks 01176 jatom_block = MODULO(tmp_block,nblocks)+1 01177 tmp_block = tmp_block/nblocks 01178 iatom_block = MODULO(tmp_block,nblocks)+1 01179 IF(jatom_block<iatom_block) CYCLE 01180 my_current_counter = my_current_counter + 1 01181 IF(my_current_counter > distribution_energy%number_of_atom_quartets) EXIT atomic_blocks 01182 01183 iatom_start = actual_x_data%blocks(iatom_block)%istart 01184 iatom_end = actual_x_data%blocks(iatom_block)%iend 01185 jatom_start = actual_x_data%blocks(jatom_block)%istart 01186 jatom_end = actual_x_data%blocks(jatom_block)%iend 01187 katom_start = actual_x_data%blocks(katom_block)%istart 01188 katom_end = actual_x_data%blocks(katom_block)%iend 01189 latom_start = actual_x_data%blocks(latom_block)%istart 01190 latom_end = actual_x_data%blocks(latom_block)%iend 01191 01192 pmax_blocks = MAX( shm_pmax_block(katom_block, iatom_block), & 01193 shm_pmax_block(latom_block, jatom_block), & 01194 shm_pmax_block(latom_block, iatom_block), & 01195 shm_pmax_block(katom_block, jatom_block) ) 01196 01197 IF( 2.0_dp*coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz ) CYCLE 01198 01199 01200 CALL build_pair_list(natom, list_ij,set_list_ij,iatom_start, iatom_end, & 01201 jatom_start, jatom_end, & 01202 kind_of,basis_parameter,particle_set, & 01203 do_periodic,screen_coeffs_set,screen_coeffs_kind, & 01204 coeffs_kind_max0,log10_eps_schwarz,cell, pmax_blocks, & 01205 shm_atomic_pair_list) 01206 01207 01208 CALL build_pair_list(natom, list_kl,set_list_kl, katom_start, katom_end, & 01209 latom_start, latom_end, & 01210 kind_of,basis_parameter,particle_set, & 01211 do_periodic,screen_coeffs_set,screen_coeffs_kind, & 01212 coeffs_kind_max0,log10_eps_schwarz, cell, pmax_blocks, & 01213 shm_atomic_pair_list) 01214 01215 DO i_list_ij=1,list_ij%n_element 01216 01217 iatom=list_ij%elements(i_list_ij)%pair(1) 01218 jatom=list_ij%elements(i_list_ij)%pair(2) 01219 i_set_list_ij_start=list_ij%elements(i_list_ij)%set_bounds(1) 01220 i_set_list_ij_stop=list_ij%elements(i_list_ij)%set_bounds(2) 01221 ikind=list_ij%elements(i_list_ij)%kind_pair(1) 01222 jkind=list_ij%elements(i_list_ij)%kind_pair(2) 01223 ra=list_ij%elements(i_list_ij)%r1 01224 rb=list_ij%elements(i_list_ij)%r2 01225 rab2=list_ij%elements(i_list_ij)%dist2 01226 01227 la_max => basis_parameter(ikind)%lmax 01228 la_min => basis_parameter(ikind)%lmin 01229 npgfa => basis_parameter(ikind)%npgf 01230 nseta = basis_parameter(ikind)%nset 01231 zeta => basis_parameter(ikind)%zet 01232 nsgfa => basis_parameter(ikind)%nsgf 01233 sphi_a_ext => basis_parameter(ikind)%sphi_ext(:,:,:,:) 01234 nsgfl_a => basis_parameter(ikind)%nsgfl 01235 sphi_a_u1=UBOUND(sphi_a_ext,1) 01236 sphi_a_u2=UBOUND(sphi_a_ext,2) 01237 sphi_a_u3=UBOUND(sphi_a_ext,3) 01238 01239 lb_max => basis_parameter(jkind)%lmax 01240 lb_min => basis_parameter(jkind)%lmin 01241 npgfb => basis_parameter(jkind)%npgf 01242 nsetb = basis_parameter(jkind)%nset 01243 zetb => basis_parameter(jkind)%zet 01244 nsgfb => basis_parameter(jkind)%nsgf 01245 sphi_b_ext => basis_parameter(jkind)%sphi_ext(:,:,:,:) 01246 nsgfl_b => basis_parameter(jkind)%nsgfl 01247 sphi_b_u1=UBOUND(sphi_b_ext,1) 01248 sphi_b_u2=UBOUND(sphi_b_ext,2) 01249 sphi_b_u3=UBOUND(sphi_b_ext,3) 01250 01251 DO i_list_kl=1,list_kl%n_element 01252 katom=list_kl%elements(i_list_kl)%pair(1) 01253 latom=list_kl%elements(i_list_kl)%pair(2) 01254 01255 IF (.NOT. (katom+latom<=iatom+jatom)) CYCLE 01256 IF( ((iatom+jatom).EQ.(katom+latom) ) .AND.(katom<iatom)) CYCLE 01257 01258 i_set_list_kl_start=list_kl%elements(i_list_kl)%set_bounds(1) 01259 i_set_list_kl_stop=list_kl%elements(i_list_kl)%set_bounds(2) 01260 kkind=list_kl%elements(i_list_kl)%kind_pair(1) 01261 lkind=list_kl%elements(i_list_kl)%kind_pair(2) 01262 rc=list_kl%elements(i_list_kl)%r1 01263 rd=list_kl%elements(i_list_kl)%r2 01264 rcd2=list_kl%elements(i_list_kl)%dist2 01265 01266 IF( do_p_screening) THEN 01267 pmax_atom = MAX(shm_pmax_atom(katom,iatom), & 01268 shm_pmax_atom(latom,jatom), & 01269 shm_pmax_atom(latom,iatom), & 01270 shm_pmax_atom(katom,jatom) ) 01271 ELSE 01272 pmax_atom = 0.0_dp 01273 END IF 01274 01275 screen_kind_ij = screen_coeffs_kind(jkind,ikind)%x(1)*rab2+& 01276 screen_coeffs_kind(jkind,ikind)%x(2) 01277 screen_kind_kl = screen_coeffs_kind(lkind,kkind)%x(1)*rcd2+& 01278 screen_coeffs_kind(lkind,kkind)%x(2) 01279 01280 IF( screen_kind_ij + screen_kind_kl + pmax_atom < log10_eps_schwarz ) CYCLE 01281 01282 !! we want to be consistent with the KS matrix. If none of the atomic indices 01283 !! is associated cycle 01284 IF(.NOT. (shm_is_assoc_atomic_block(latom,iatom)>=1 .AND. & 01285 shm_is_assoc_atomic_block(katom,iatom)>=1 .AND. & 01286 shm_is_assoc_atomic_block(katom,jatom)>=1 .AND. & 01287 shm_is_assoc_atomic_block(latom,jatom)>=1 ) ) CYCLE 01288 01289 !! calculate symmetry_factor accordin to degeneracy of atomic quartet 01290 symm_fac = 0.5_dp 01291 IF(iatom==jatom) symm_fac = symm_fac*2.0_dp 01292 IF(katom==latom) symm_fac = symm_fac*2.0_dp 01293 IF(iatom==katom .AND. jatom==latom .AND. iatom/=jatom .AND. katom/=latom) symm_fac = symm_fac*2.0_dp 01294 IF(iatom==katom .AND. iatom==jatom .AND. katom==latom) symm_fac = symm_fac*2.0_dp 01295 symm_fac = 1.0_dp / symm_fac 01296 01297 lc_max => basis_parameter(kkind)%lmax 01298 lc_min => basis_parameter(kkind)%lmin 01299 npgfc => basis_parameter(kkind)%npgf 01300 zetc => basis_parameter(kkind)%zet 01301 nsgfc => basis_parameter(kkind)%nsgf 01302 sphi_c_ext => basis_parameter(kkind)%sphi_ext(:,:,:,:) 01303 nsgfl_c => basis_parameter(kkind)%nsgfl 01304 sphi_c_u1=UBOUND(sphi_c_ext,1) 01305 sphi_c_u2=UBOUND(sphi_c_ext,2) 01306 sphi_c_u3=UBOUND(sphi_c_ext,3) 01307 01308 ld_max => basis_parameter(lkind)%lmax 01309 ld_min => basis_parameter(lkind)%lmin 01310 npgfd => basis_parameter(lkind)%npgf 01311 zetd => basis_parameter(lkind)%zet 01312 nsgfd => basis_parameter(lkind)%nsgf 01313 sphi_d_ext => basis_parameter(lkind)%sphi_ext(:,:,:,:) 01314 nsgfl_d => basis_parameter(lkind)%nsgfl 01315 sphi_d_u1=UBOUND(sphi_d_ext,1) 01316 sphi_d_u2=UBOUND(sphi_d_ext,2) 01317 sphi_d_u3=UBOUND(sphi_d_ext,3) 01318 01319 atomic_offset_bd = shm_atomic_block_offset(jatom,latom) 01320 atomic_offset_bc = shm_atomic_block_offset(jatom,katom) 01321 atomic_offset_ad = shm_atomic_block_offset(iatom,latom) 01322 atomic_offset_ac = shm_atomic_block_offset(iatom,katom) 01323 01324 IF( jatom < latom) THEN 01325 offset_bd_set => shm_set_offset(:,:,lkind,jkind) 01326 ELSE 01327 offset_bd_set => shm_set_offset(:,:,jkind,lkind) 01328 END IF 01329 IF( jatom < katom) THEN 01330 offset_bc_set => shm_set_offset(:,:,kkind,jkind) 01331 ELSE 01332 offset_bc_set => shm_set_offset(:,:,jkind,kkind) 01333 END IF 01334 IF( iatom < latom) THEN 01335 offset_ad_set => shm_set_offset(:,:,lkind,ikind) 01336 ELSE 01337 offset_ad_set => shm_set_offset(:,:,ikind,lkind) 01338 END IF 01339 IF( iatom < katom) THEN 01340 offset_ac_set => shm_set_offset(:,:,kkind,ikind) 01341 ELSE 01342 offset_ac_set => shm_set_offset(:,:,ikind,kkind) 01343 END IF 01344 01345 IF( do_p_screening) THEN 01346 swap_id = 0 01347 kind_kind_idx = INT(get_1D_idx(kkind,ikind,INT(nkind,int_8))) 01348 IF( ikind >= kkind) THEN 01349 ptr_p_1 => shm_initial_p(kind_kind_idx)%p_kind(:,:, & 01350 actual_x_data%map_atom_to_kind_atom(katom), & 01351 actual_x_data%map_atom_to_kind_atom(iatom)) 01352 ELSE 01353 ptr_p_1 => shm_initial_p(kind_kind_idx)%p_kind(:,:, & 01354 actual_x_data%map_atom_to_kind_atom(iatom), & 01355 actual_x_data%map_atom_to_kind_atom(katom)) 01356 swap_id = swap_id + 1 01357 END IF 01358 kind_kind_idx = INT(get_1D_idx(lkind,jkind,INT(nkind,int_8))) 01359 IF( jkind>=lkind ) THEN 01360 ptr_p_2 => shm_initial_p(kind_kind_idx)%p_kind(:,:, & 01361 actual_x_data%map_atom_to_kind_atom(latom), & 01362 actual_x_data%map_atom_to_kind_atom(jatom)) 01363 ELSE 01364 ptr_p_2 => shm_initial_p(kind_kind_idx)%p_kind(:,:, & 01365 actual_x_data%map_atom_to_kind_atom(jatom), & 01366 actual_x_data%map_atom_to_kind_atom(latom)) 01367 swap_id = swap_id + 2 01368 END IF 01369 kind_kind_idx = INT(get_1D_idx(lkind,ikind,INT(nkind,int_8))) 01370 IF( ikind>=lkind ) THEN 01371 ptr_p_3 => shm_initial_p(kind_kind_idx)%p_kind(:,:, & 01372 actual_x_data%map_atom_to_kind_atom(latom), & 01373 actual_x_data%map_atom_to_kind_atom(iatom)) 01374 ELSE 01375 ptr_p_3 => shm_initial_p(kind_kind_idx)%p_kind(:,:, & 01376 actual_x_data%map_atom_to_kind_atom(iatom), & 01377 actual_x_data%map_atom_to_kind_atom(latom)) 01378 swap_id = swap_id + 4 01379 END IF 01380 kind_kind_idx = INT(get_1D_idx(kkind,jkind,INT(nkind,int_8))) 01381 IF( jkind>=kkind) THEN 01382 ptr_p_4 => shm_initial_p(kind_kind_idx)%p_kind(:,:, & 01383 actual_x_data%map_atom_to_kind_atom(katom), & 01384 actual_x_data%map_atom_to_kind_atom(jatom)) 01385 ELSE 01386 ptr_p_4 => shm_initial_p(kind_kind_idx)%p_kind(:,:, & 01387 actual_x_data%map_atom_to_kind_atom(jatom), & 01388 actual_x_data%map_atom_to_kind_atom(katom)) 01389 swap_id = swap_id + 8 01390 END IF 01391 END IF 01392 01393 01394 !! At this stage, check for memory used in compression 01395 01396 IF(geometry_did_change ) THEN 01397 IF(.NOT.memory_parameter%do_all_on_the_fly) THEN 01398 ! ** We know the maximum amount of integrals that we can store per MPI-process 01399 ! ** Now we need to sum the current memory usage among all openMP threads 01400 ! ** We can just read what is currently stored on the corresponding x_data type 01401 ! ** If this is thread i and it tries to read the data from thread j, that is 01402 ! ** currently writing that data, we just dont care, because the possible error 01403 ! ** is of the order of CACHE_SIZE 01404 mem_compression_counter = 0 01405 DO i=1,n_threads 01406 mem_compression_counter = mem_compression_counter + & 01407 x_data(irep,i)%memory_parameter%actual_memory_usage *& 01408 memory_parameter%cache_size 01409 END DO 01410 IF(mem_compression_counter > memory_parameter%max_compression_counter) THEN 01411 buffer_overflow = .TRUE. 01412 IF( do_dynamic_load_balancing ) THEN 01413 distribution_energy%ram_counter = counter 01414 ELSE 01415 memory_parameter%ram_counter = counter 01416 END IF 01417 ELSE 01418 counter = counter + 1 01419 buffer_overflow = .FALSE. 01420 END IF 01421 END IF 01422 ELSE 01423 IF(.NOT.memory_parameter%do_all_on_the_fly) THEN 01424 IF( do_dynamic_load_balancing) THEN 01425 IF(distribution_energy%ram_counter == counter) THEN 01426 buffer_overflow = .TRUE. 01427 ELSE 01428 counter = counter + 1 01429 buffer_overflow = .FALSE. 01430 END IF 01431 01432 ELSE 01433 IF(memory_parameter%ram_counter == counter) THEN 01434 buffer_overflow = .TRUE. 01435 ELSE 01436 counter = counter + 1 01437 buffer_overflow = .FALSE. 01438 END IF 01439 END IF 01440 END IF 01441 END IF 01442 01443 IF( buffer_overflow .AND. do_disk_storage ) THEN 01444 use_disk_storage = .TRUE. 01445 buffer_overflow = .FALSE. 01446 END IF 01447 01448 IF( use_disk_storage ) THEN 01449 mem_compression_counter_disk = memory_parameter%actual_memory_usage_disk *& 01450 memory_parameter%cache_size 01451 IF(mem_compression_counter_disk > memory_parameter%max_compression_counter_disk) THEN 01452 buffer_overflow = .TRUE. 01453 use_disk_storage = .FALSE. 01454 END IF 01455 END IF 01456 01457 DO i_set_list_ij=i_set_list_ij_start, i_set_list_ij_stop 01458 iset=set_list_ij(i_set_list_ij)%pair(1) 01459 jset=set_list_ij(i_set_list_ij)%pair(2) 01460 01461 ncob = npgfb(jset)*ncoset(lb_max(jset)) 01462 max_val1 = screen_coeffs_set(jset,iset,jkind,ikind)%x(1)*rab2 + & 01463 screen_coeffs_set(jset,iset,jkind,ikind)%x(2) 01464 01465 IF( max_val1 + screen_kind_kl + pmax_atom < log10_eps_schwarz) CYCLE 01466 01467 sphi_a_ext_set => sphi_a_ext(:,:,:,iset) 01468 sphi_b_ext_set => sphi_b_ext(:,:,:,jset) 01469 DO i_set_list_kl=i_set_list_kl_start, i_set_list_kl_stop 01470 kset=set_list_kl(i_set_list_kl)%pair(1) 01471 lset=set_list_kl(i_set_list_kl)%pair(2) 01472 01473 max_val2_set = (screen_coeffs_set(lset,kset,lkind,kkind)%x(1)*rcd2 + & 01474 screen_coeffs_set(lset,kset,lkind,kkind)%x(2) ) 01475 max_val2 = max_val1 + max_val2_set 01476 01477 !! Near field screening 01478 IF(max_val2 + pmax_atom <log10_eps_schwarz) CYCLE 01479 sphi_c_ext_set => sphi_c_ext(:,:,:,kset) 01480 sphi_d_ext_set => sphi_d_ext(:,:,:,lset) 01481 !! get max_vals if we screen on initial density 01482 IF( do_p_screening) THEN 01483 CALL get_pmax_val(ptr_p_1, ptr_p_2, ptr_p_3, ptr_p_4, & 01484 iset, jset, kset, lset, & 01485 pmax_entry, swap_id) 01486 ELSE 01487 pmax_entry = 0.0_dp 01488 END IF 01489 log10_pmax = pmax_entry 01490 max_val2 = max_val2 + log10_pmax 01491 IF(max_val2<log10_eps_schwarz) CYCLE 01492 pmax_entry = EXP(log10_pmax*ln_10) 01493 01494 !! store current number of integrals, update total number and number of integrals in buffer 01495 current_counter = nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset) 01496 IF(buffer_overflow) THEN 01497 neris_onthefly = neris_onthefly + current_counter 01498 END IF 01499 01500 !! Get integrals from buffer and update Kohn-Sham matrix 01501 IF(.NOT.buffer_overflow .AND. .NOT.geometry_did_change) THEN 01502 nints = current_counter 01503 IF(.NOT. use_disk_storage) THEN 01504 CALL hfx_get_single_cache_element(estimate_to_store_int, 6,& 01505 maxval_cache, maxval_container, memory_parameter%actual_memory_usage, & 01506 use_disk_storage) 01507 ELSE 01508 CALL hfx_get_single_cache_element(estimate_to_store_int, 6,& 01509 maxval_cache_disk, maxval_container_disk, & 01510 memory_parameter%actual_memory_usage_disk, & 01511 use_disk_storage) 01512 END IF 01513 spherical_estimate = SET_EXPONENT(1.0_dp,estimate_to_store_int+1) 01514 IF(spherical_estimate*pmax_entry<eps_schwarz) CYCLE 01515 nbits = EXPONENT(ANINT(spherical_estimate*pmax_entry/eps_storage)) + 1 01516 buffer_left = nints 01517 buffer_start = 1 01518 DO WHILE (buffer_left > 0 ) 01519 buffer_size = MIN(buffer_left, cache_size) 01520 IF(.NOT. use_disk_storage) THEN 01521 CALL hfx_get_mult_cache_elements(primitive_integrals(buffer_start), & 01522 buffer_size, nbits, & 01523 integral_caches(nbits), & 01524 integral_containers(nbits), & 01525 eps_storage, pmax_entry, & 01526 memory_parameter%actual_memory_usage, & 01527 use_disk_storage) 01528 ELSE 01529 CALL hfx_get_mult_cache_elements(primitive_integrals(buffer_start), & 01530 buffer_size, nbits, & 01531 integral_caches_disk(nbits), & 01532 integral_containers_disk(nbits), & 01533 eps_storage, pmax_entry, & 01534 memory_parameter%actual_memory_usage_disk, & 01535 use_disk_storage) 01536 END IF 01537 buffer_left = buffer_left - buffer_size 01538 buffer_start = buffer_start + buffer_size 01539 END DO 01540 END IF 01541 !! Calculate integrals if we run out of buffer or the geometry did change 01542 IF(geometry_did_change .OR. buffer_overflow) THEN 01543 max_contraction_val = max_contraction(iset,iatom) * & 01544 max_contraction(jset,jatom) * & 01545 max_contraction(kset,katom) * & 01546 max_contraction(lset,latom) * pmax_entry 01547 tmp_R_1 => radii_pgf(:,:,jset,iset,jkind,ikind) 01548 tmp_R_2 => radii_pgf(:,:,lset,kset,lkind,kkind) 01549 tmp_screen_pgf1 => screen_coeffs_pgf(:,:,jset,iset,jkind,ikind) 01550 tmp_screen_pgf2 => screen_coeffs_pgf(:,:,lset,kset,lkind,kkind) 01551 01552 CALL coulomb4(private_lib, ra, rb, rc, rd, npgfa(iset), npgfb(jset), npgfc(kset), npgfd(lset), & 01553 la_min(iset), la_max(iset), lb_min(jset), lb_max(jset),& 01554 lc_min(kset), lc_max(kset), ld_min(lset), ld_max(lset),& 01555 nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset),& 01556 sphi_a_u1,sphi_a_u2,sphi_a_u3,& 01557 sphi_b_u1,sphi_b_u2,sphi_b_u3,& 01558 sphi_c_u1,sphi_c_u2,sphi_c_u3,& 01559 sphi_d_u1,sphi_d_u2,sphi_d_u3,& 01560 zeta(1:npgfa(iset),iset), zetb(1:npgfb(jset),jset),& 01561 zetc(1:npgfc(kset),kset), zetd(1:npgfd(lset),lset),& 01562 primitive_integrals,& 01563 potential_parameter, & 01564 actual_x_data%neighbor_cells, screen_coeffs_set(jset,iset,jkind,ikind)%x,& 01565 screen_coeffs_set(lset,kset,lkind,kkind)%x, eps_schwarz, & 01566 max_contraction_val, cartesian_estimate, cell, neris_tmp,& 01567 log10_pmax, log10_eps_schwarz, & 01568 tmp_R_1, tmp_R_2, tmp_screen_pgf1, tmp_screen_pgf2,& 01569 pgf_list_ij, pgf_list_kl, pgf_product_list, & 01570 nsgfl_a(:,iset), nsgfl_b(:,jset),& 01571 nsgfl_c(:,kset), nsgfl_d(:,lset),& 01572 sphi_a_ext_set,& 01573 sphi_b_ext_set,& 01574 sphi_c_ext_set,& 01575 sphi_d_ext_set,& 01576 ee_work,ee_work2,ee_buffer1,ee_buffer2, ee_primitives_tmp,& 01577 nimages,do_periodic, p_work) 01578 01579 nints = nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset) 01580 neris_total = neris_total + nints 01581 nprim_ints = nprim_ints + neris_tmp 01582 IF(cartesian_estimate == 0.0_dp) cartesian_estimate = TINY(cartesian_estimate) 01583 estimate_to_store_int = EXPONENT(cartesian_estimate) 01584 estimate_to_store_int = MAX(estimate_to_store_int,-15_int_8) 01585 cartesian_estimate = SET_EXPONENT(1.0_dp,estimate_to_store_int+1) 01586 IF(.NOT.buffer_overflow .AND. geometry_did_change) THEN 01587 IF(cartesian_estimate < eps_schwarz ) THEN 01588 IF(.NOT. use_disk_storage) THEN 01589 CALL hfx_add_single_cache_element(estimate_to_store_int, 6,& 01590 maxval_cache, maxval_container, memory_parameter%actual_memory_usage, & 01591 use_disk_storage, max_val_memory, error) 01592 ELSE 01593 CALL hfx_add_single_cache_element(estimate_to_store_int, 6,& 01594 maxval_cache_disk, maxval_container_disk, & 01595 memory_parameter%actual_memory_usage_disk, & 01596 use_disk_storage, error=error) 01597 END IF 01598 END IF 01599 END IF 01600 01601 IF(cartesian_estimate<eps_schwarz) CYCLE 01602 01603 !! Compress the array for storage 01604 spherical_estimate = 0.0_dp 01605 DO i=1,nints 01606 spherical_estimate = MAX(spherical_estimate,ABS(primitive_integrals(i))) 01607 END DO 01608 01609 IF(spherical_estimate == 0.0_dp) spherical_estimate = TINY(spherical_estimate) 01610 estimate_to_store_int = EXPONENT(spherical_estimate) 01611 estimate_to_store_int = MAX(estimate_to_store_int, -15_int_8) 01612 01613 IF(.NOT. buffer_overflow .AND. geometry_did_change) THEN 01614 IF(.NOT.use_disk_storage) THEN 01615 CALL hfx_add_single_cache_element(estimate_to_store_int, 6, & 01616 maxval_cache, maxval_container, & 01617 memory_parameter%actual_memory_usage, & 01618 use_disk_storage, max_val_memory, error) 01619 ELSE 01620 CALL hfx_add_single_cache_element(estimate_to_store_int, 6, & 01621 maxval_cache_disk, maxval_container_disk, & 01622 memory_parameter%actual_memory_usage_disk, & 01623 use_disk_storage, error=error) 01624 END IF 01625 END IF 01626 spherical_estimate = SET_EXPONENT(1.0_dp,estimate_to_store_int+1) 01627 IF(spherical_estimate*pmax_entry<eps_schwarz) CYCLE 01628 IF(.NOT.buffer_overflow) THEN 01629 nbits = EXPONENT(ANINT(spherical_estimate*pmax_entry/eps_storage)) + 1 01630 buffer_left = nints 01631 buffer_start = 1 01632 IF( .NOT. use_disk_storage) THEN 01633 neris_incore = neris_incore + nints 01634 ELSE 01635 neris_disk = neris_disk + nints 01636 END IF 01637 DO WHILE (buffer_left > 0 ) 01638 buffer_size = MIN(buffer_left, CACHE_SIZE) 01639 IF(.NOT. use_disk_storage) THEN 01640 CALL hfx_add_mult_cache_elements(primitive_integrals(buffer_start),& 01641 buffer_size, nbits, & 01642 integral_caches(nbits), & 01643 integral_containers(nbits), & 01644 eps_storage, pmax_entry, & 01645 memory_parameter%actual_memory_usage, & 01646 use_disk_storage, error) 01647 ELSE 01648 CALL hfx_add_mult_cache_elements(primitive_integrals(buffer_start),& 01649 buffer_size, nbits, & 01650 integral_caches_disk(nbits), & 01651 integral_containers_disk(nbits), & 01652 eps_storage, pmax_entry, & 01653 memory_parameter%actual_memory_usage_disk, & 01654 use_disk_storage, error) 01655 END IF 01656 buffer_left = buffer_left - buffer_size 01657 buffer_start = buffer_start + buffer_size 01658 END DO 01659 ELSE 01660 !! In order to be consistent with in-core part, round all the eris wrt. eps_schwarz 01661 DO i=1,nints 01662 primitive_integrals(i) = primitive_integrals(i)*pmax_entry 01663 IF(ABS(primitive_integrals(i)) > eps_storage) THEN 01664 primitive_integrals(i) = ANINT(primitive_integrals(i)/eps_storage, dp)*eps_storage/pmax_entry 01665 ELSE 01666 primitive_integrals(i) = 0.0_dp 01667 END IF 01668 END DO 01669 END IF 01670 END IF 01671 !! Update Kohn-Sham matrix 01672 CALL update_fock_matrix(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), & 01673 fac, symm_fac, full_density, full_ks, & 01674 primitive_integrals, pbd_buf, pbc_buf, pad_buf, pac_buf, kbd_buf,& 01675 kbc_buf, kad_buf, kac_buf, iatom, jatom, katom, latom,& 01676 iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set, offset_ac_set,& 01677 atomic_offset_bd, atomic_offset_bc, atomic_offset_ad, atomic_offset_ac) 01678 IF( .NOT. treat_lsd_in_core ) THEN 01679 IF( nspins == 2 ) THEN 01680 CALL update_fock_matrix(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), & 01681 fac, symm_fac, full_density_beta, full_ks_beta, & 01682 primitive_integrals, pbd_buf, pbc_buf, pad_buf, pac_buf, kbd_buf,& 01683 kbc_buf, kad_buf, kac_buf, iatom, jatom, katom, latom,& 01684 iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set,offset_ac_set,& 01685 atomic_offset_bd, atomic_offset_bc, atomic_offset_ad, atomic_offset_ac) 01686 END IF 01687 END IF 01688 END DO ! i_set_list_kl 01689 END DO ! i_set_list_ij 01690 IF( do_disk_storage ) THEN 01691 buffer_overflow = .TRUE. 01692 END IF 01693 END DO !i_list_ij 01694 END DO !i_list_kl 01695 END DO atomic_blocks 01696 bintime_stop = m_walltime() 01697 IF (geometry_did_change) THEN 01698 distribution_energy%time_first_scf=bintime_stop-bintime_start 01699 ELSE 01700 distribution_energy%time_other_scf= & 01701 distribution_energy%time_other_scf+bintime_stop-bintime_start 01702 ENDIF 01703 !$OMP MASTER 01704 CALL timestop(handle_bin) 01705 !$OMP END MASTER 01706 END DO !bin 01707 01708 01709 !$OMP MASTER 01710 logger => cp_error_get_logger(error) 01711 do_print_load_balance_info = .FALSE. 01712 do_print_load_balance_info = BTEST(cp_print_key_should_output(logger%iter_info,hfx_section,& 01713 "LOAD_BALANCE%PRINT/LOAD_BALANCE_INFO",error=error),cp_p_file) 01714 !$OMP END MASTER 01715 !$OMP BARRIER 01716 IF(do_print_load_balance_info) THEN 01717 iw = -1 01718 !$OMP MASTER 01719 iw = cp_print_key_unit_nr(logger,hfx_section,"LOAD_BALANCE%PRINT/LOAD_BALANCE_INFO",& 01720 extension=".scfLog",error=error) 01721 !$OMP END MASTER 01722 01723 CALL collect_load_balance_info(para_env, actual_x_data, iw, n_threads, i_thread, & 01724 hfx_do_eval_energy, error) 01725 01726 !$OMP MASTER 01727 CALL cp_print_key_finished_output(iw,logger,hfx_section,& 01728 "LOAD_BALANCE%PRINT/LOAD_BALANCE_INFO", error=error) 01729 !$OMP END MASTER 01730 END IF 01731 01732 !$OMP BARRIER 01733 !$OMP MASTER 01734 memsize_after=m_memory() 01735 !$OMP END MASTER 01736 !$OMP BARRIER 01737 01738 DEALLOCATE(primitive_integrals,STAT=stat) 01739 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01740 !$OMP BARRIER 01741 !! Get some number about ERIS 01742 !$OMP ATOMIC 01743 shm_neris_total = shm_neris_total + neris_total 01744 !$OMP ATOMIC 01745 shm_neris_onthefly = shm_neris_onthefly + neris_onthefly 01746 !$OMP ATOMIC 01747 shm_nprim_ints = shm_nprim_ints + nprim_ints 01748 01749 storage_counter_integrals = memory_parameter%actual_memory_usage * & 01750 memory_parameter%cache_size 01751 stor_count_int_disk = memory_parameter%actual_memory_usage_disk * & 01752 memory_parameter%cache_size 01753 stor_count_max_val = max_val_memory * memory_parameter%cache_size 01754 !$OMP ATOMIC 01755 shm_storage_counter_integrals = shm_storage_counter_integrals + storage_counter_integrals 01756 !$OMP ATOMIC 01757 shm_stor_count_int_disk = shm_stor_count_int_disk + stor_count_int_disk 01758 !$OMP ATOMIC 01759 shm_neris_incore = shm_neris_incore + neris_incore 01760 !$OMP ATOMIC 01761 shm_neris_disk = shm_neris_disk + neris_disk 01762 !$OMP ATOMIC 01763 shm_stor_count_max_val = shm_stor_count_max_val + stor_count_max_val 01764 !$OMP BARRIER 01765 01766 01767 ! ** Calculate how much memory has already been used (might be needed for in-core forces 01768 !$OMP MASTER 01769 shm_mem_compression_counter = 0 01770 DO i=1,n_threads 01771 shm_mem_compression_counter = shm_mem_compression_counter + & 01772 x_data(irep,i)%memory_parameter%actual_memory_usage *& 01773 memory_parameter%cache_size 01774 END DO 01775 !$OMP END MASTER 01776 !$OMP BARRIER 01777 actual_x_data%memory_parameter%final_comp_counter_energy = shm_mem_compression_counter 01778 01779 01780 !$OMP MASTER 01781 !! Calculate the exchange energies from the Kohn-Sham matrix. Before we can go on, we have to symmetrize. 01782 ene_x_aa = 0.0_dp 01783 ene_x_bb = 0.0_dp 01784 01785 mb_size_p = shm_block_offset(ncpu+1) /1024/128 01786 mb_size_f = shm_block_offset(ncpu+1) /1024/128 01787 IF (.NOT. treat_lsd_in_core ) THEN 01788 IF( nspins == 2) THEN 01789 mb_size_f = mb_size_f * 2 01790 mb_size_p = mb_size_p * 2 01791 END IF 01792 END IF 01793 !! size of primitive_integrals(not shared) 01794 mb_size_buffers = INT(nsgf_max,int_8)**4 * n_threads 01795 !! fock density buffers (not shared) 01796 mb_size_buffers = mb_size_buffers + INT(nsgf_max,int_8)**2 * n_threads 01797 subtr_size_mb = subtr_size_mb + 8_int_8*nsgf_max**2*n_threads 01798 !! size of screening functions (shared) 01799 mb_size_buffers = mb_size_buffers + max_pgf**2*max_set**2*nkind**2 & 01800 + max_set**2*nkind**2 & 01801 + nkind**2 & 01802 + max_pgf**2*max_set**2*nkind**2 01803 !! is_assoc (shared) 01804 mb_size_buffers = mb_size_buffers+natom**2 01805 ! ** pmax_atom (shared) 01806 IF( do_p_screening) THEN 01807 mb_size_buffers = mb_size_buffers + natom**2 01808 END IF 01809 IF(screening_parameter%do_p_screening_forces) THEN 01810 IF( memory_parameter%treat_forces_in_core) THEN 01811 mb_size_buffers = mb_size_buffers + natom**2 01812 END IF 01813 END IF 01814 ! ** Initial P only MAX(alpha,beta) (shared) 01815 IF(do_p_screening .OR. screening_parameter%do_p_screening_forces) THEN 01816 mb_size_buffers = mb_size_buffers + memory_parameter%size_p_screen 01817 END IF 01818 ! ** In core forces require their own initial P 01819 IF(screening_parameter%do_p_screening_forces) THEN 01820 IF( memory_parameter%treat_forces_in_core) THEN 01821 mb_size_buffers = mb_size_buffers + memory_parameter%size_p_screen 01822 END IF 01823 END IF 01824 01825 !! mb 01826 mb_size_buffers = mb_size_buffers /1024/128 01827 01828 CALL timestop(handle_main) 01829 01830 ene_x_aa_diag = 0.0_dp 01831 ene_x_bb_diag = 0.0_dp 01832 DO iatom = 1,natom 01833 ikind = kind_of(iatom) 01834 nseta = basis_parameter(ikind)%nset 01835 nsgfa => basis_parameter(ikind)%nsgf 01836 jatom = iatom 01837 jkind = kind_of(jatom) 01838 nsetb = basis_parameter(jkind)%nset 01839 nsgfb => basis_parameter(jkind)%nsgf 01840 act_atomic_block_offset = shm_atomic_block_offset(jatom,iatom) 01841 DO iset=1,nseta 01842 DO jset=1,nsetb 01843 act_set_offset = shm_set_offset(jset,iset,jkind,ikind) 01844 i = act_set_offset + act_atomic_block_offset - 1 01845 DO ma=1,nsgfa(iset) 01846 j = shm_set_offset(iset,jset,jkind,ikind) + act_atomic_block_offset - 1 + ma -1 01847 DO mb=1,nsgfb(jset) 01848 IF(i>j) THEN 01849 full_ks(i) = (full_ks(i) + full_ks(j)) 01850 full_ks(j) = full_ks(i) 01851 IF( .NOT. treat_lsd_in_core .AND. nspins == 2 ) THEN 01852 full_ks_beta(i) = (full_ks_beta(i) + full_ks_beta(j)) 01853 full_ks_beta(j) = full_ks_beta(i) 01854 END IF 01855 END IF 01856 ! ** For adiabatically rescaled functionals we need the energy coming from the diagonal elements 01857 IF( i==j ) THEN 01858 ene_x_aa_diag = ene_x_aa_diag + full_ks(i)*full_density(i) 01859 IF( .NOT. treat_lsd_in_core .AND. nspins == 2 ) THEN 01860 ene_x_bb_diag = ene_x_bb_diag + full_ks_beta(i)*full_density_beta(i) 01861 END IF 01862 END IF 01863 i = i+ 1 01864 j = j+nsgfa(iset) 01865 END DO 01866 END DO 01867 END DO 01868 END DO 01869 END DO 01870 01871 CALL mp_sync(para_env%group) 01872 IF( distribute_fock_matrix) THEN 01873 !! Distribute the current KS-matrix to all the processes 01874 CALL timeset(routineN//"_dist_KS",handle_dist_ks) 01875 CALL distribute_ks_matrix(para_env, full_ks, ks_matrix(ispin)%matrix, shm_number_of_p_entries, & 01876 shm_block_offset, natom, kind_of, basis_parameter, & 01877 off_diag_fac=0.5_dp,error=error) 01878 01879 NULLIFY(full_ks) 01880 DEALLOCATE(shm_master_x_data%full_ks_alpha, STAT=stat) 01881 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01882 IF( .NOT. treat_lsd_in_core ) THEN 01883 IF(nspins == 2) THEN 01884 CALL distribute_ks_matrix(para_env, full_ks_beta, ks_matrix(2)%matrix, shm_number_of_p_entries, & 01885 shm_block_offset, natom, kind_of, basis_parameter, & 01886 off_diag_fac=0.5_dp, error=error) 01887 NULLIFY(full_ks_beta) 01888 DEALLOCATE(shm_master_x_data%full_ks_beta, STAT=stat) 01889 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01890 END IF 01891 END IF 01892 CALL timestop(handle_dist_ks) 01893 END IF 01894 01895 01896 IF( distribute_fock_matrix) THEN 01897 !! ** Calculate the exchange energy 01898 CALL cp_dbcsr_trace(ks_matrix(ispin)%matrix, rho_ao(ispin)%matrix, & 01899 ene_x_aa, error=error) 01900 IF( nspins == 2 .AND. .NOT. treat_lsd_in_core ) THEN 01901 CALL cp_dbcsr_trace(ks_matrix(2)%matrix, rho_ao(2)%matrix, & 01902 ene_x_bb, error=error) 01903 END IF 01904 01905 !! Update energy type 01906 energy%ex = 0.5_dp*(ene_x_aa+ene_x_bb) + energy%ex 01907 ELSE 01908 ! ** It is easier to correct the following expression by the diagonal energy contribution, 01909 ! ** than explicitly going throuhg the diagonal elements 01910 DO pa=1,SIZE(full_ks) 01911 ene_x_aa = ene_x_aa + full_ks(pa)*full_density(pa) 01912 END DO 01913 ! ** Now correct 01914 ene_x_aa = (ene_x_aa + ene_x_aa_diag)*0.5_dp 01915 IF( nspins == 2 ) THEN 01916 DO pa=1,SIZE(full_ks) 01917 ene_x_bb = ene_x_bb + full_ks_beta(pa)*full_density_beta(pa) 01918 END DO 01919 ! ** Now correct 01920 ene_x_bb = (ene_x_bb + ene_x_bb_diag)*0.5_dp 01921 END IF 01922 CALL mp_sum(ene_x_aa, para_env%group) 01923 IF(nspins==2) CALL mp_sum(ene_x_bb,para_env%group) 01924 energy%ex = 0.5_dp*(ene_x_aa+ene_x_bb) 01925 END IF 01926 01927 !! Print some memeory information if this is the first step 01928 IF( geometry_did_change ) THEN 01929 tmp_i8(1:8)=(/shm_storage_counter_integrals,shm_neris_onthefly,shm_neris_incore,shm_neris_disk, & 01930 shm_neris_total,shm_stor_count_int_disk,shm_nprim_ints,shm_stor_count_max_val/) 01931 CALL mp_sum(tmp_i8,para_env%group) 01932 shm_storage_counter_integrals=tmp_i8(1) 01933 shm_neris_onthefly=tmp_i8(2) 01934 shm_neris_incore=tmp_i8(3) 01935 shm_neris_disk=tmp_i8(4) 01936 shm_neris_total=tmp_i8(5) 01937 shm_stor_count_int_disk=tmp_i8(6) 01938 shm_nprim_ints=tmp_i8(7) 01939 shm_stor_count_max_val=tmp_i8(8) 01940 CALL mp_max(memsize_after, para_env%group) 01941 mem_eris = (shm_storage_counter_integrals+128*1024-1)/1024/128 01942 compression_factor = REAL(shm_neris_incore,dp)/REAL(shm_storage_counter_integrals,dp) 01943 mem_eris_disk = (shm_stor_count_int_disk+128*1024-1)/1024/128 01944 compression_factor_disk = REAL(shm_neris_disk,dp)/REAL(shm_stor_count_int_disk,dp) 01945 mem_max_val = (shm_stor_count_max_val+128*1024-1)/1024/128 01946 01947 IF( shm_neris_incore == 0 ) THEN 01948 mem_eris = 0 01949 compression_factor = 0.0_dp 01950 END IF 01951 IF( shm_neris_disk == 0 ) THEN 01952 mem_eris_disk = 0 01953 compression_factor_disk = 0.0_dp 01954 END IF 01955 01956 iw = cp_print_key_unit_nr(logger,hfx_section,"HF_INFO",& 01957 extension=".scfLog",error=error) 01958 IF (iw>0) THEN 01959 WRITE (UNIT=iw,FMT="((T3,A,T60,I21))")& 01960 "HFX_MEM_INFO| Number of cart. primitive ERI's calculated: ",shm_nprim_ints 01961 01962 WRITE (UNIT=iw,FMT="((T3,A,T60,I21))")& 01963 "HFX_MEM_INFO| Number of sph. ERI's calculated: ",shm_neris_total 01964 01965 WRITE (UNIT=iw,FMT="((T3,A,T60,I21))")& 01966 "HFX_MEM_INFO| Number of sph. ERI's stored in-core: ",shm_neris_incore 01967 01968 WRITE (UNIT=iw,FMT="((T3,A,T60,I21))")& 01969 "HFX_MEM_INFO| Number of sph. ERI's stored on disk: ",shm_neris_disk 01970 01971 WRITE (UNIT=iw,FMT="((T3,A,T60,I21))")& 01972 "HFX_MEM_INFO| Number of sph. ERI's calculated on the fly: ",shm_neris_onthefly 01973 01974 WRITE (UNIT=iw,FMT="((T3,A,T60,I21))")& 01975 "HFX_MEM_INFO| Total memory consumption ERI's RAM [MB's]: ",mem_eris 01976 01977 WRITE (UNIT=iw,FMT="((T3,A,T60,I21))")& 01978 "HFX_MEM_INFO| Whereof max-vals [MB's]: ",mem_max_val 01979 01980 WRITE (UNIT=iw,FMT="((T3,A,T60,F21.2))")& 01981 "HFX_MEM_INFO| Total compression factor ERI's RAM: ",compression_factor 01982 01983 WRITE (UNIT=iw,FMT="((T3,A,T60,I21))")& 01984 "HFX_MEM_INFO| Total memory consumption ERI's disk [MB's]: ",mem_eris_disk 01985 01986 WRITE (UNIT=iw,FMT="((T3,A,T60,F21.2))")& 01987 "HFX_MEM_INFO| Total compression factor ERI's disk: ",compression_factor_disk 01988 01989 WRITE (UNIT=iw,FMT="((T3,A,T60,I21))")& 01990 "HFX_MEM_INFO| Size of density/Fock matrix [MB's]: ", 2_int_8*mb_size_p 01991 01992 IF( do_periodic) THEN 01993 WRITE (UNIT=iw,FMT="((T3,A,T60,I21))")& 01994 "HFX_MEM_INFO| Size of buffers [MB's]: ", mb_size_buffers 01995 WRITE (UNIT=iw,FMT="((T3,A,T60,I21))")& 01996 "HFX_MEM_INFO| Number of periodic image cells considered: ", SIZE(shm_master_x_data%neighbor_cells) 01997 ELSE 01998 WRITE (UNIT=iw,FMT="((T3,A,T60,I21))")& 01999 "HFX_MEM_INFO| Size of buffers [MB's]: ", mb_size_buffers 02000 END IF 02001 WRITE (UNIT=iw,FMT="((T3,A,T60,I21),/)")& 02002 "HFX_MEM_INFO| Est. max. program size after HFX [MB's]:", memsize_after/(1024*1024) 02003 02004 END IF 02005 02006 CALL cp_print_key_finished_output(iw,logger,hfx_section,& 02007 "HF_INFO", error=error) 02008 END IF 02009 !$OMP END MASTER 02010 02011 !! flush caches if the geometry changed 02012 IF( do_dynamic_load_balancing ) THEN 02013 my_bin_size = SIZE(actual_x_data%distribution_energy) 02014 ELSE 02015 my_bin_size = 1 02016 END IF 02017 02018 IF(geometry_did_change ) THEN 02019 IF(.NOT.memory_parameter%do_all_on_the_fly) THEN 02020 DO bin = 1,my_bin_size 02021 maxval_cache => actual_x_data%maxval_cache(bin) 02022 maxval_container => actual_x_data%maxval_container(bin) 02023 integral_caches => actual_x_data%integral_caches(:,bin) 02024 integral_containers => actual_x_data%integral_containers(:,bin) 02025 CALL hfx_flush_last_cache(bits_max_val, maxval_cache, maxval_container, memory_parameter%actual_memory_usage, & 02026 .FALSE., error) 02027 DO i=1,64 02028 CALL hfx_flush_last_cache(i, integral_caches(i), integral_containers(i), & 02029 memory_parameter%actual_memory_usage, .FALSE., error) 02030 END DO 02031 END DO 02032 END IF 02033 END IF 02034 !! reset all caches except we calculate all on the fly 02035 IF(.NOT.memory_parameter%do_all_on_the_fly) THEN 02036 DO bin = 1,my_bin_size 02037 maxval_cache => actual_x_data%maxval_cache(bin) 02038 maxval_container => actual_x_data%maxval_container(bin) 02039 integral_caches => actual_x_data%integral_caches(:,bin) 02040 integral_containers => actual_x_data%integral_containers(:,bin) 02041 02042 CALL hfx_reset_cache_and_container(maxval_cache, maxval_container, memory_parameter%actual_memory_usage, .FALSE.) 02043 DO i=1,64 02044 CALL hfx_reset_cache_and_container(integral_caches(i), integral_containers(i), & 02045 memory_parameter%actual_memory_usage, & 02046 .FALSE.) 02047 END DO 02048 END DO 02049 END IF 02050 02051 !! Since the I/O routines are no thread-safe, i.e. the procedure to get the unit number, put a lock here 02052 !$OMP CRITICAL 02053 IF(do_disk_storage ) THEN 02054 !! flush caches if the geometry changed 02055 IF(geometry_did_change ) THEN 02056 CALL hfx_flush_last_cache(bits_max_val, maxval_cache_disk, maxval_container_disk, & 02057 memory_parameter%actual_memory_usage_disk, .TRUE., error) 02058 DO i=1,64 02059 CALL hfx_flush_last_cache(i, integral_caches_disk(i), integral_containers_disk(i), & 02060 memory_parameter%actual_memory_usage_disk, .TRUE., error) 02061 END DO 02062 END IF 02063 !! reset all caches except we calculate all on the fly 02064 CALL hfx_reset_cache_and_container(maxval_cache_disk, maxval_container_disk, memory_parameter%actual_memory_usage_disk, & 02065 do_disk_storage) 02066 DO i=1,64 02067 CALL hfx_reset_cache_and_container(integral_caches_disk(i), integral_containers_disk(i),& 02068 memory_parameter%actual_memory_usage_disk, do_disk_storage) 02069 END DO 02070 END IF 02071 !$OMP END CRITICAL 02072 !$OMP BARRIER 02073 !! Clean up 02074 DEALLOCATE(last_sgf_global,STAT=stat) 02075 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02076 !$OMP MASTER 02077 DEALLOCATE(full_density,STAT=stat) 02078 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02079 IF( .NOT. treat_lsd_in_core ) THEN 02080 IF(nspins==2) THEN 02081 DEALLOCATE(full_density_beta,STAT=stat) 02082 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02083 END IF 02084 END IF 02085 IF (do_dynamic_load_balancing) THEN 02086 DEALLOCATE(shm_master_x_data%task_list, STAT=stat) 02087 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02088 END IF 02089 !$OMP END MASTER 02090 DEALLOCATE(pbd_buf, pbc_buf, pad_buf, pac_buf, STAT=stat) 02091 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02092 DEALLOCATE(kbd_buf,kbc_buf,kad_buf,kac_buf, STAT=stat) 02093 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02094 DEALLOCATE(set_list_ij, set_list_kl, STAT=stat) 02095 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02096 02097 DO i=1,max_pgf**2 02098 DEALLOCATE(pgf_list_ij(i)%image_list, STAT=stat) 02099 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02100 DEALLOCATE(pgf_list_kl(i)%image_list, STAT=stat) 02101 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02102 END DO 02103 02104 DEALLOCATE(pgf_list_ij, STAT=stat) 02105 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02106 DEALLOCATE(pgf_list_kl, STAT=stat) 02107 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02108 DEALLOCATE(pgf_product_list, STAT=stat) 02109 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02110 02111 DEALLOCATE(max_contraction, kind_of, STAT=stat) 02112 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02113 02114 DEALLOCATE(ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp, STAT=stat) 02115 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02116 02117 IF( .NOT. has_iso_c_binding) THEN 02118 DEALLOCATE(p_work, STAT=stat) 02119 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02120 END IF 02121 02122 DEALLOCATE(nimages, STAT=stat) 02123 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02124 02125 !$OMP BARRIER 02126 !$OMP END PARALLEL 02127 02128 #if defined(__USE_PAT) 02129 CALL PAT_sampling_state(PAT_STATE_OFF,istat) 02130 #endif 02131 02132 CALL timestop(handle) 02133 END SUBROUTINE integrate_four_center 02134 02135 ! ***************************************************************************** 02164 SUBROUTINE coulomb4(lib, ra, rb, rc, rd, npgfa, npgfb, npgfc, npgfd, & 02165 la_min, la_max, lb_min, lb_max,& 02166 lc_min, lc_max, ld_min, ld_max, nsgfa, nsgfb, nsgfc, nsgfd,& 02167 sphi_a_u1,sphi_a_u2,sphi_a_u3,& 02168 sphi_b_u1,sphi_b_u2,sphi_b_u3,& 02169 sphi_c_u1,sphi_c_u2,sphi_c_u3,& 02170 sphi_d_u1,sphi_d_u2,sphi_d_u3,& 02171 zeta, zetb, zetc, zetd,& 02172 primitive_integrals,& 02173 potential_parameter, neighbor_cells,& 02174 screen1, screen2, eps_schwarz, max_contraction_val,& 02175 cart_estimate, cell, neris_tmp, log10_pmax,& 02176 log10_eps_schwarz,R1_pgf, R2_pgf, pgf1, pgf2,& 02177 pgf_list_ij, pgf_list_kl,& 02178 pgf_product_list,& 02179 nsgfl_a, nsgfl_b, nsgfl_c,& 02180 nsgfl_d,& 02181 sphi_a_ext, sphi_b_ext, sphi_c_ext, sphi_d_ext,& 02182 ee_work,ee_work2,ee_buffer1,ee_buffer2, ee_primitives_tmp,& 02183 nimages, do_periodic, p_work) 02184 02185 TYPE(lib_int) :: lib 02186 REAL(dp), INTENT(IN) :: ra(3), rb(3), rc(3), rd(3) 02187 INTEGER, INTENT(IN) :: npgfa, npgfb, npgfc, npgfd, la_min, la_max, 02188 lb_min, lb_max, lc_min, lc_max, ld_min, ld_max, nsgfa, nsgfb, nsgfc, 02189 nsgfd, sphi_a_u1, sphi_a_u2, sphi_a_u3, sphi_b_u1, sphi_b_u2, 02190 sphi_b_u3, sphi_c_u1, sphi_c_u2, sphi_c_u3, sphi_d_u1, sphi_d_u2, 02191 sphi_d_u3 02192 REAL(dp), DIMENSION(1:npgfa), INTENT(IN) :: zeta 02193 REAL(dp), DIMENSION(1:npgfb), INTENT(IN) :: zetb 02194 REAL(dp), DIMENSION(1:npgfc), INTENT(IN) :: zetc 02195 REAL(dp), DIMENSION(1:npgfd), INTENT(IN) :: zetd 02196 REAL(dp), 02197 DIMENSION(nsgfa, nsgfb, nsgfc, nsgfd) :: primitive_integrals 02198 TYPE(hfx_potential_type) :: potential_parameter 02199 TYPE(hfx_cell_type), DIMENSION(:), 02200 POINTER :: neighbor_cells 02201 REAL(dp), INTENT(IN) :: screen1(2), screen2(2), 02202 eps_schwarz, 02203 max_contraction_val 02204 REAL(dp) :: cart_estimate 02205 TYPE(cell_type), POINTER :: cell 02206 INTEGER(int_8) :: neris_tmp 02207 REAL(dp), INTENT(IN) :: log10_pmax, log10_eps_schwarz 02208 TYPE(hfx_screen_coeff_type), 02209 DIMENSION(:, :), POINTER :: R1_pgf, R2_pgf, pgf1, pgf2 02210 TYPE(hfx_pgf_list), DIMENSION(*) :: pgf_list_ij, pgf_list_kl 02211 TYPE(hfx_pgf_product_list), DIMENSION(*) :: pgf_product_list 02212 INTEGER, DIMENSION(0:), INTENT(IN) :: nsgfl_a, nsgfl_b, nsgfl_c, 02213 nsgfl_d 02214 REAL(dp), INTENT(IN) :: sphi_a_ext(sphi_a_u1,sphi_a_u2,sphi_a_u3), 02215 sphi_b_ext(sphi_b_u1,sphi_b_u2,sphi_b_u3), 02216 sphi_c_ext(sphi_c_u1,sphi_c_u2,sphi_c_u3), 02217 sphi_d_ext(sphi_d_u1,sphi_d_u2,sphi_d_u3) 02218 REAL(dp), DIMENSION(*) :: ee_work, ee_work2, 02219 ee_buffer1, ee_buffer2, 02220 ee_primitives_tmp 02221 INTEGER, DIMENSION(*) :: nimages 02222 LOGICAL, INTENT(IN) :: do_periodic 02223 REAL(dp), DIMENSION(:), POINTER :: p_work 02224 02225 INTEGER :: ipgf, jpgf, kpgf, la, lb, lc, ld, list_ij, list_kl, lpgf, 02226 max_l, ncoa, ncob, ncoc, ncod, nelements_ij, nelements_kl, nneighbors3, 02227 nproducts, nsgfla, nsgflb, nsgflc, nsgfld, nsoa, nsob, nsoc, nsod, 02228 s_offset_a, s_offset_b, s_offset_c, s_offset_d 02229 REAL(dp) :: EtaInv, tmp_max, ZetaInv 02230 02231 CALL build_pair_list_pgf(npgfa, npgfb, pgf_list_ij, zeta, zetb, screen1, screen2,& 02232 pgf1, R1_pgf, log10_pmax, log10_eps_schwarz, ra, rb,& 02233 nelements_ij, & 02234 neighbor_cells, nimages, do_periodic) 02235 CALL build_pair_list_pgf(npgfc, npgfd, pgf_list_kl, zetc, zetd, screen2, screen1, & 02236 pgf2, R2_pgf, log10_pmax, log10_eps_schwarz, rc, rd,& 02237 nelements_kl, & 02238 neighbor_cells, nimages, do_periodic) 02239 02240 cart_estimate = 0.0_dp 02241 neris_tmp = 0 02242 primitive_integrals = 0.0_dp 02243 max_l = la_max + lb_max + lc_max + ld_max 02244 DO list_ij = 1,nelements_ij 02245 ZetaInv = pgf_list_ij(list_ij)%ZetaInv 02246 ipgf = pgf_list_ij(list_ij)%ipgf 02247 jpgf = pgf_list_ij(list_ij)%jpgf 02248 02249 DO list_kl = 1,nelements_kl 02250 EtaInv = pgf_list_kl(list_kl)%ZetaInv 02251 kpgf = pgf_list_kl(list_kl)%ipgf 02252 lpgf = pgf_list_kl(list_kl)%jpgf 02253 02254 nneighbors3 = SIZE(neighbor_cells)*pgf_list_ij(list_ij)%nimages*pgf_list_kl(list_kl)%nimages 02255 CALL build_pgf_product_list(nneighbors3, pgf_list_ij(list_ij), pgf_list_kl(list_kl), pgf_product_list, & 02256 nproducts, log10_pmax, log10_eps_schwarz, neighbor_cells, cell, & 02257 potential_parameter, max_l, do_periodic) 02258 02259 s_offset_a = 0 02260 DO la = la_min,la_max 02261 s_offset_b = 0 02262 ncoa = nco(la) 02263 nsgfla = nsgfl_a(la) 02264 nsoa = nso(la) 02265 02266 DO lb = lb_min, lb_max 02267 s_offset_c = 0 02268 ncob = nco(lb) 02269 nsgflb = nsgfl_b(lb) 02270 nsob = nso(lb) 02271 02272 DO lc = lc_min, lc_max 02273 s_offset_d = 0 02274 ncoc = nco(lc) 02275 nsgflc = nsgfl_c(lc) 02276 nsoc = nso(lc) 02277 02278 DO ld = ld_min, ld_max 02279 ncod = nco(ld) 02280 nsgfld = nsgfl_d(ld) 02281 nsod = nso(ld) 02282 02283 tmp_max = 0.0_dp 02284 CALL evaluate_eri(lib, nproducts, pgf_product_list, & 02285 la, lb, lc ,ld,& 02286 ncoa,ncob,ncoc,ncod,& 02287 nsgfa, nsgfb, nsgfc, nsgfd,& 02288 primitive_integrals, & 02289 max_contraction_val, tmp_max, eps_schwarz,& 02290 neris_tmp,ZetaInv,EtaInv,& 02291 s_offset_a, s_offset_b, s_offset_c, s_offset_d,& 02292 nsgfla, nsgflb, nsgflc, nsgfld, nsoa, nsob, nsoc, nsod, & 02293 sphi_a_ext(1,la+1,ipgf), & 02294 sphi_b_ext(1,lb+1,jpgf), & 02295 sphi_c_ext(1,lc+1,kpgf), & 02296 sphi_d_ext(1,ld+1,lpgf), & 02297 ee_work,ee_work2,ee_buffer1,ee_buffer2, ee_primitives_tmp,& 02298 p_work) 02299 cart_estimate = MAX(tmp_max,cart_estimate) 02300 s_offset_d = s_offset_d + nsod*nsgfld 02301 END DO !ld 02302 s_offset_c = s_offset_c + nsoc*nsgflc 02303 END DO !lc 02304 s_offset_b = s_offset_b + nsob*nsgflb 02305 END DO !lb 02306 s_offset_a = s_offset_a + nsoa*nsgfla 02307 END DO !la 02308 END DO 02309 END DO 02310 02311 END SUBROUTINE coulomb4 02312 02313 ! ***************************************************************************** 02324 PURE FUNCTION get_1D_idx(i,j,N) 02325 INTEGER, INTENT(IN) :: i, j 02326 INTEGER(int_8), INTENT(IN) :: N 02327 INTEGER(int_8) :: get_1D_idx 02328 02329 INTEGER(int_8) :: min_ij 02330 02331 min_ij = MIN(i,j) 02332 get_1D_idx = min_ij*N + MAX(i,j) - (min_ij-1)*min_ij/2 - N 02333 02334 END FUNCTION get_1D_idx 02335 02336 ! ***************************************************************************** 02353 02354 SUBROUTINE update_fock_matrix(ma_max, mb_max, mc_max, md_max, & 02355 fac, symm_fac, density, ks, prim, & 02356 pbd, pbc, pad, pac, kbd, kbc, kad, kac, & 02357 iatom, jatom, katom, latom, & 02358 iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set,& 02359 offset_ac_set, atomic_offset_bd, atomic_offset_bc, atomic_offset_ad,& 02360 atomic_offset_ac) 02361 02362 INTEGER, INTENT(IN) :: ma_max, mb_max, mc_max, md_max 02363 REAL(dp), INTENT(IN) :: fac, symm_fac 02364 REAL(dp), DIMENSION(:), INTENT(IN) :: density 02365 REAL(dp), DIMENSION(:), INTENT(INOUT) :: ks 02366 REAL(dp), DIMENSION(ma_max*mb_max*mc_max& *md_max), INTENT(IN) :: prim 02367 REAL(dp), DIMENSION(*), INTENT(INOUT) :: pbd, pbc, pad, pac, kbd, kbc, 02368 kad, kac 02369 INTEGER, INTENT(IN) :: iatom, jatom, katom, latom, 02370 iset, jset, kset, lset 02371 INTEGER, DIMENSION(:, :), POINTER :: offset_bd_set, offset_bc_set, 02372 offset_ad_set, offset_ac_set 02373 INTEGER, INTENT(IN) :: atomic_offset_bd, 02374 atomic_offset_bc, 02375 atomic_offset_ad, 02376 atomic_offset_ac 02377 02378 INTEGER :: i, j, ma, mb, mc, md, 02379 offset_ac, offset_ad, 02380 offset_bc, offset_bd 02381 02382 IF( jatom >= latom ) THEN 02383 i = 1 02384 offset_bd = offset_bd_set(jset,lset) + atomic_offset_bd - 1 02385 j = offset_bd 02386 DO md=1,md_max 02387 DO mb=1,mb_max 02388 pbd(i) = density(j) 02389 i = i + 1 02390 j = j + 1 02391 END DO 02392 END DO 02393 ELSE 02394 i = 1 02395 offset_bd = offset_bd_set(lset,jset) + atomic_offset_bd - 1 02396 DO md=1,md_max 02397 j = offset_bd + md - 1 02398 DO mb=1,mb_max 02399 pbd(i) = density(j) 02400 i = i + 1 02401 j = j + md_max 02402 END DO 02403 END DO 02404 END IF 02405 IF ( jatom >= katom ) THEN 02406 i = 1 02407 offset_bc = offset_bc_set(jset,kset) + atomic_offset_bc - 1 02408 j = offset_bc 02409 DO mc=1,mc_max 02410 DO mb=1,mb_max 02411 pbc(i) = density(j) 02412 i = i + 1 02413 j = j + 1 02414 END DO 02415 END DO 02416 ELSE 02417 i = 1 02418 offset_bc = offset_bc_set(kset,jset) + atomic_offset_bc - 1 02419 DO mc=1,mc_max 02420 j = offset_bc + mc -1 02421 DO mb=1,mb_max 02422 pbc(i) = density(j) 02423 i = i + 1 02424 j = j + mc_max 02425 END DO 02426 END DO 02427 END IF 02428 IF( iatom >= latom) THEN 02429 i = 1 02430 offset_ad = offset_ad_set(iset,lset) + atomic_offset_ad - 1 02431 j = offset_ad 02432 DO md=1,md_max 02433 DO ma=1,ma_max 02434 pad(i) = density(j) 02435 i = i + 1 02436 j = j + 1 02437 END DO 02438 END DO 02439 ELSE 02440 i = 1 02441 offset_ad = offset_ad_set(lset,iset) + atomic_offset_ad - 1 02442 DO md=1,md_max 02443 j = offset_ad + md -1 02444 DO ma=1,ma_max 02445 pad(i) = density(j) 02446 i = i + 1 02447 j = j + md_max 02448 END DO 02449 END DO 02450 END IF 02451 IF( iatom >= katom ) THEN 02452 i = 1 02453 offset_ac = offset_ac_set(iset,kset) + atomic_offset_ac - 1 02454 j = offset_ac 02455 DO mc=1,mc_max 02456 DO ma=1,ma_max 02457 pac(i) = density(j) 02458 i = i + 1 02459 j = j + 1 02460 END DO 02461 END DO 02462 ELSE 02463 i = 1 02464 offset_ac = offset_ac_set(kset,iset) + atomic_offset_ac - 1 02465 DO mc=1,mc_max 02466 j = offset_ac + mc - 1 02467 DO ma=1,ma_max 02468 pac(i) = density(j) 02469 i = i + 1 02470 j = j + mc_max 02471 END DO 02472 END DO 02473 END IF 02474 02475 CALL contract_block(ma_max,mb_max,mc_max,md_max,kbd,kbc,kad,kac,pbd,pbc,pad,pac,prim,fac*symm_fac) 02476 02477 IF(jatom>=latom) THEN 02478 i = 1 02479 j = offset_bd 02480 DO md=1,md_max 02481 DO mb=1,mb_max 02482 !$OMP ATOMIC 02483 ks(j) = ks(j) + kbd(i) 02484 i = i + 1 02485 j = j + 1 02486 END DO 02487 END DO 02488 ELSE 02489 i = 1 02490 DO md=1,md_max 02491 j = offset_bd+md-1 02492 DO mb=1,mb_max 02493 !$OMP ATOMIC 02494 ks(j) = ks(j) + kbd(i) 02495 i = i + 1 02496 j = j + md_max 02497 END DO 02498 END DO 02499 END IF 02500 IF(jatom>=katom) THEN 02501 i = 1 02502 j = offset_bc 02503 DO mc=1,mc_max 02504 DO mb=1,mb_max 02505 !$OMP ATOMIC 02506 ks(j) = ks(j) + kbc(i) 02507 i = i + 1 02508 j = j + 1 02509 END DO 02510 END DO 02511 ELSE 02512 i = 1 02513 DO mc=1,mc_max 02514 j = offset_bc + mc -1 02515 DO mb=1,mb_max 02516 !$OMP ATOMIC 02517 ks(j) = ks(j) + kbc(i) 02518 i = i + 1 02519 j = j + mc_max 02520 END DO 02521 END DO 02522 END IF 02523 IF(iatom>=latom) THEN 02524 i = 1 02525 j = offset_ad 02526 DO md=1,md_max 02527 DO ma=1,ma_max 02528 !$OMP ATOMIC 02529 ks(j) = ks(j) + kad(i) 02530 i = i + 1 02531 j = j + 1 02532 END DO 02533 END DO 02534 ELSE 02535 i = 1 02536 DO md=1,md_max 02537 j = offset_ad + md - 1 02538 DO ma=1,ma_max 02539 !$OMP ATOMIC 02540 ks(j) = ks(j) + kad(i) 02541 i = i + 1 02542 j = j + md_max 02543 END DO 02544 END DO 02545 END IF 02546 IF(iatom>=katom) THEN 02547 i = 1 02548 j = offset_ac 02549 DO mc=1,mc_max 02550 DO ma=1,ma_max 02551 !$OMP ATOMIC 02552 ks(j) = ks(j) + kac(i) 02553 i = i + 1 02554 j = j + 1 02555 END DO 02556 END DO 02557 ELSE 02558 i = 1 02559 DO mc=1,mc_max 02560 j = offset_ac + mc -1 02561 DO ma=1,ma_max 02562 !$OMP ATOMIC 02563 ks(j) = ks(j) + kac(i) 02564 i = i + 1 02565 j = j + mc_max 02566 END DO 02567 END DO 02568 END IF 02569 END SUBROUTINE update_fock_matrix 02570 02571 #include "hfx_get_pmax_val.f90" 02572 END MODULE hfx_energy_potential 02573
1.7.3