|
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 admm_methods 00013 USE admm_types, ONLY: admm_create_block_list,& 00014 admm_env_create,& 00015 admm_type 00016 USE cp_control_types, ONLY: admm_control_type 00017 USE cp_dbcsr_interface, ONLY: & 00018 cp_dbcsr_add, cp_dbcsr_col_block_sizes, cp_dbcsr_copy, & 00019 cp_dbcsr_create, cp_dbcsr_distribution, cp_dbcsr_get_block_p, & 00020 cp_dbcsr_get_data_size, cp_dbcsr_get_data_type, & 00021 cp_dbcsr_get_num_blocks, cp_dbcsr_init, cp_dbcsr_iterator_blocks_left, & 00022 cp_dbcsr_iterator_next_block, cp_dbcsr_iterator_start, & 00023 cp_dbcsr_iterator_stop, cp_dbcsr_row_block_sizes, cp_dbcsr_scale, & 00024 cp_dbcsr_set 00025 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 00026 copy_fm_to_dbcsr,& 00027 cp_dbcsr_alloc_block_from_nbl,& 00028 cp_dbcsr_deallocate_matrix,& 00029 cp_dbcsr_plus_fm_fm_t 00030 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix 00031 USE cp_dbcsr_types, ONLY: cp_dbcsr_iterator,& 00032 cp_dbcsr_p_type,& 00033 cp_dbcsr_type 00034 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,& 00035 cp_fm_gemm,& 00036 cp_fm_scale,& 00037 cp_fm_scale_and_add,& 00038 cp_fm_schur_product,& 00039 cp_fm_upper_to_full 00040 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,& 00041 cp_fm_cholesky_invert,& 00042 cp_fm_cholesky_reduce,& 00043 cp_fm_cholesky_restore 00044 USE cp_fm_diag, ONLY: cp_fm_syevd 00045 USE cp_fm_types, ONLY: cp_fm_p_type,& 00046 cp_fm_set_all,& 00047 cp_fm_set_element,& 00048 cp_fm_to_fm,& 00049 cp_fm_type 00050 USE cp_output_handling, ONLY: cp_p_file,& 00051 cp_print_key_finished_output,& 00052 cp_print_key_should_output,& 00053 cp_print_key_unit_nr 00054 USE cp_para_types, ONLY: cp_para_env_type 00055 USE dbcsr_types, ONLY: dbcsr_type_no_symmetry,& 00056 dbcsr_type_symmetric 00057 USE input_constants, ONLY: & 00058 do_admm_basis_set_projection, do_admm_block_aux_basis_off, & 00059 do_admm_block_aux_basis_on, do_admm_block_density_matrix, & 00060 do_admm_block_purify_blocked, do_admm_block_purify_full, & 00061 do_admm_block_purify_off, do_admm_purify_cauchy, & 00062 do_admm_purify_cauchy_subspace, do_admm_purify_mo_diag, & 00063 do_admm_purify_mo_no_diag, do_admm_purify_none, & 00064 do_hfx_potential_coulomb, do_hfx_potential_short, & 00065 do_hfx_potential_truncated, use_aux_fit_basis_set, use_orb_basis_set, & 00066 xc_funct_no_shortcut 00067 USE input_section_types, ONLY: section_vals_duplicate,& 00068 section_vals_get_subs_vals,& 00069 section_vals_get_subs_vals2,& 00070 section_vals_remove_values,& 00071 section_vals_type,& 00072 section_vals_val_get,& 00073 section_vals_val_set 00074 USE kinds, ONLY: dp 00075 USE mathconstants 00076 USE particle_types, ONLY: particle_type 00077 USE qs_environment_types, ONLY: get_qs_env,& 00078 qs_environment_type 00079 USE qs_mo_methods, ONLY: calculate_density_matrix 00080 USE qs_mo_types, ONLY: get_mo_set,& 00081 mo_set_p_type,& 00082 mo_set_type 00083 USE qs_overlap, ONLY: build_overlap_matrix 00084 USE timings, ONLY: timeset,& 00085 timestop 00086 #include "cp_common_uses.h" 00087 00088 IMPLICIT NONE 00089 PRIVATE 00090 00091 PUBLIC admm_fit_mo_coeffs, & 00092 admm_merge_ks_matrix, & 00093 admm_merge_mo_derivs, & 00094 remove_ks_matrix, & 00095 calc_aux_mo_derivs_none, & 00096 calc_mixed_overlap_force, & 00097 create_admm_xc_section,& 00098 print_matlab_matrix,& 00099 admm_calculate_density_matrix,& 00100 fit_mo_coeffs_diag 00101 00102 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'admm_methods' 00103 00104 !*** 00105 00106 CONTAINS 00107 00108 SUBROUTINE admm_calculate_density_matrix(admm_env,mo_set,density_matrix, density_matrix_aux,& 00109 ispin, nspins, use_dbcsr,error) 00110 TYPE(admm_type), POINTER :: admm_env 00111 TYPE(mo_set_type), POINTER :: mo_set 00112 TYPE(cp_dbcsr_type), POINTER :: density_matrix, 00113 density_matrix_aux 00114 INTEGER :: ispin, nspins 00115 LOGICAL, INTENT(IN), OPTIONAL :: use_dbcsr 00116 TYPE(cp_error_type), INTENT(inout) :: error 00117 00118 CHARACTER(len=*), PARAMETER :: 00119 routineN = 'admm_calculate_density_matrix', 00120 routineP = moduleN//':'//routineN 00121 00122 INTEGER :: handle 00123 LOGICAL :: failure 00124 00125 failure = .FALSE. 00126 00127 CALL timeset(routineN,handle) 00128 00129 SELECT CASE(admm_env%method_id) 00130 CASE (do_admm_basis_set_projection) 00131 SELECT CASE(admm_env%purification_method) 00132 CASE(do_admm_purify_none) 00133 CALL calculate_density_matrix(mo_set, density_matrix_aux, use_dbcsr, error=error) 00134 CASE(do_admm_purify_cauchy) 00135 CALL calculate_dm_mo_no_diag(admm_env, mo_set, density_matrix_aux, ispin, error) 00136 CALL purify_density_matrix_cauchy(admm_env, mo_set, density_matrix_aux, ispin, error) 00137 CASE(do_admm_purify_cauchy_subspace) 00138 CALL calculate_dm_mo_no_diag(admm_env, mo_set, density_matrix_aux, ispin, error) 00139 CASE(do_admm_purify_mo_diag) 00140 CALL calculate_density_matrix(mo_set, density_matrix_aux, use_dbcsr, error=error) 00141 CASE(do_admm_purify_mo_no_diag) 00142 CALL calculate_dm_mo_no_diag(admm_env, mo_set, density_matrix_aux, ispin, error) 00143 END SELECT 00144 00145 CASE(do_admm_block_density_matrix) 00146 SELECT CASE(admm_env%block_purification_method) 00147 CASE(do_admm_block_purify_full) 00148 SELECT CASE(admm_env%block_projection_method) 00149 CASE (do_admm_block_aux_basis_off) 00150 CALL blockify_density_matrix(admm_env, density_matrix, density_matrix_aux, ispin, nspins, error) 00151 CALL purify_dm_cauchy_blocked(admm_env, mo_set, density_matrix_aux, ispin, error) 00152 CASE (do_admm_block_aux_basis_on) 00153 END SELECT 00154 CASE(do_admm_block_purify_off) 00155 SELECT CASE(admm_env%block_projection_method) 00156 CASE (do_admm_block_aux_basis_off) 00157 CALL blockify_density_matrix(admm_env, density_matrix, density_matrix_aux, ispin, nspins, error) 00158 CASE (do_admm_block_aux_basis_on) 00159 END SELECT 00160 CASE(do_admm_block_purify_blocked) 00161 SELECT CASE(admm_env%block_projection_method) 00162 CASE (do_admm_block_aux_basis_off) 00163 CALL blockify_density_matrix(admm_env, density_matrix, density_matrix_aux, ispin, nspins, error) 00164 CALL purify_dm_cauchy_blocked(admm_env, mo_set, density_matrix_aux, ispin, error) 00165 CASE (do_admm_block_aux_basis_on) 00166 END SELECT 00167 END SELECT 00168 END SELECT 00169 00170 00171 CALL timestop(handle) 00172 00173 END SUBROUTINE admm_calculate_density_matrix 00174 00175 SUBROUTINE admm_merge_mo_derivs(ispin, admm_env, mo_set, mo_coeff, mo_coeff_aux_fit, mo_derivs, & 00176 mo_derivs_aux_fit, matrix_ks_aux_fit, error) 00177 INTEGER, INTENT(IN) :: ispin 00178 TYPE(admm_type), POINTER :: admm_env 00179 TYPE(mo_set_type), POINTER :: mo_set 00180 TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit 00181 TYPE(cp_fm_p_type), DIMENSION(:), 00182 POINTER :: mo_derivs, mo_derivs_aux_fit 00183 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00184 POINTER :: matrix_ks_aux_fit 00185 TYPE(cp_error_type), INTENT(INOUT) :: error 00186 00187 CHARACTER(LEN=*), PARAMETER :: routineN = 'admm_merge_mo_derivs', 00188 routineP = moduleN//':'//routineN 00189 00190 INTEGER :: handle 00191 LOGICAL :: failure 00192 00193 failure = .FALSE. 00194 00195 CALL timeset(routineN,handle) 00196 00197 SELECT CASE(admm_env%method_id) 00198 CASE (do_admm_basis_set_projection) 00199 SELECT CASE(admm_env%purification_method) 00200 CASE(do_admm_purify_none) 00201 CASE(do_admm_purify_cauchy) 00202 CASE(do_admm_purify_cauchy_subspace) 00203 CASE(do_admm_purify_mo_diag) 00204 CALL merge_mo_derivs_diag(ispin, admm_env, mo_set, mo_coeff, mo_coeff_aux_fit,& 00205 mo_derivs,mo_derivs_aux_fit, matrix_ks_aux_fit,& 00206 error) 00207 CASE(do_admm_purify_mo_no_diag) 00208 CALL merge_mo_derivs_no_diag(ispin, admm_env, mo_set, mo_coeff, mo_coeff_aux_fit,& 00209 mo_derivs,mo_derivs_aux_fit, matrix_ks_aux_fit,& 00210 error) 00211 END SELECT 00212 00213 CASE(do_admm_block_density_matrix) 00214 00215 END SELECT 00216 00217 CALL timestop(handle) 00218 00219 END SUBROUTINE admm_merge_mo_derivs 00220 00221 SUBROUTINE admm_merge_ks_matrix(ispin, admm_env, mo_set, mo_coeff, mo_coeff_aux_fit, & 00222 matrix_ks, matrix_ks_aux_fit, matrix_s, & 00223 matrix_p_aux_fit, matrix_p, error) 00224 INTEGER, INTENT(IN) :: ispin 00225 TYPE(admm_type), POINTER :: admm_env 00226 TYPE(mo_set_type), POINTER :: mo_set 00227 TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit 00228 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00229 POINTER :: matrix_ks, matrix_ks_aux_fit, 00230 matrix_s, matrix_p_aux_fit, 00231 matrix_p 00232 TYPE(cp_error_type), INTENT(INOUT) :: error 00233 00234 CHARACTER(LEN=*), PARAMETER :: routineN = 'admm_merge_ks_matrix', 00235 routineP = moduleN//':'//routineN 00236 00237 INTEGER :: blk, handle, iatom, jatom 00238 LOGICAL :: failure 00239 REAL(dp), DIMENSION(:, :), POINTER :: sparse_block 00240 TYPE(cp_dbcsr_iterator) :: iter 00241 00242 failure = .FALSE. 00243 00244 CALL timeset(routineN,handle) 00245 00246 SELECT CASE(admm_env%method_id) 00247 CASE (do_admm_basis_set_projection) 00248 SELECT CASE(admm_env%purification_method) 00249 CASE(do_admm_purify_none) 00250 CALL merge_ks_matrix_none(ispin, admm_env, mo_set, mo_coeff,& 00251 mo_coeff_aux_fit, matrix_ks, matrix_ks_aux_fit, & 00252 matrix_s, matrix_p_aux_fit, error) 00253 CASE(do_admm_purify_cauchy) 00254 CALL merge_ks_matrix_cauchy(ispin, admm_env, mo_set, mo_coeff,& 00255 mo_coeff_aux_fit, matrix_ks, matrix_ks_aux_fit, & 00256 matrix_s, matrix_p_aux_fit, matrix_p, error) 00257 CASE(do_admm_purify_cauchy_subspace) 00258 CALL merge_ks_matrix_cauchy_subspace(ispin, admm_env, mo_set, mo_coeff,& 00259 mo_coeff_aux_fit, matrix_ks, matrix_ks_aux_fit, & 00260 matrix_s, matrix_p_aux_fit, error) 00261 CASE(do_admm_purify_mo_diag) 00262 CASE(do_admm_purify_mo_no_diag) 00263 END SELECT 00264 00265 CASE(do_admm_block_density_matrix) 00266 SELECT CASE(admm_env%block_purification_method) 00267 CASE(do_admm_block_purify_full) 00268 SELECT CASE(admm_env%block_projection_method) 00269 CASE (do_admm_block_aux_basis_off) 00270 CALL merge_ks_matrix_cauchy_blocked(ispin, admm_env, mo_set, mo_coeff,& 00271 mo_coeff_aux_fit, matrix_ks, matrix_ks_aux_fit, & 00272 matrix_s, matrix_p_aux_fit, matrix_p, error) 00273 CASE (do_admm_block_aux_basis_on) 00274 END SELECT 00275 CASE(do_admm_block_purify_off) 00276 SELECT CASE(admm_env%block_projection_method) 00277 CASE (do_admm_block_aux_basis_off) 00278 CALL cp_dbcsr_iterator_start(iter, matrix_ks_aux_fit(ispin)%matrix) 00279 DO WHILE (cp_dbcsr_iterator_blocks_left (iter)) 00280 CALL cp_dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block, blk) 00281 IF( admm_env%block_map(iatom,jatom) == 0 ) THEN 00282 sparse_block = 0.0_dp 00283 END IF 00284 END DO 00285 CALL cp_dbcsr_iterator_stop(iter) 00286 CALL cp_dbcsr_add(matrix_ks(ispin)%matrix, matrix_ks_aux_fit(ispin)%matrix, 1.0_dp, 1.0_dp, error) 00287 CASE (do_admm_block_aux_basis_on) 00288 END SELECT 00289 CASE(do_admm_block_purify_blocked) 00290 SELECT CASE(admm_env%block_projection_method) 00291 CASE (do_admm_block_aux_basis_off) 00292 CALL merge_ks_matrix_cauchy_blocked(ispin, admm_env, mo_set, mo_coeff,& 00293 mo_coeff_aux_fit, matrix_ks, matrix_ks_aux_fit, & 00294 matrix_s, matrix_p_aux_fit, matrix_p, error) 00295 CASE (do_admm_block_aux_basis_on) 00296 END SELECT 00297 END SELECT 00298 END SELECT 00299 00300 CALL timestop(handle) 00301 00302 END SUBROUTINE admm_merge_ks_matrix 00303 00304 SUBROUTINE admm_fit_mo_coeffs(qs_env, admm_env, admm_control, para_env, matrix_s_aux_fit, matrix_s_mixed, & 00305 mos, mos_aux_fit, geometry_did_change, error) 00306 00307 TYPE(qs_environment_type), POINTER :: qs_env 00308 TYPE(admm_type), POINTER :: admm_env 00309 TYPE(admm_control_type), POINTER :: admm_control 00310 TYPE(cp_para_env_type), POINTER :: para_env 00311 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00312 POINTER :: matrix_s_aux_fit, 00313 matrix_s_mixed 00314 TYPE(mo_set_p_type), DIMENSION(:), 00315 POINTER :: mos, mos_aux_fit 00316 LOGICAL, INTENT(IN) :: geometry_did_change 00317 TYPE(cp_error_type), INTENT(inout) :: error 00318 00319 CHARACTER(LEN=*), PARAMETER :: routineN = 'admm_fit_mo_coeffs', 00320 routineP = moduleN//':'//routineN 00321 00322 INTEGER :: handle, natom 00323 TYPE(particle_type), DIMENSION(:), 00324 POINTER :: particle_set 00325 TYPE(section_vals_type), POINTER :: admm_block_section, input, 00326 xc_section 00327 00328 CALL timeset(routineN,handle) 00329 00330 NULLIFY(xc_section, admm_block_section) 00331 00332 IF (.NOT.(ASSOCIATED(admm_env) )) THEN 00333 CALL admm_env_create(mos, mos_aux_fit, & 00334 para_env, admm_env,& 00335 error) 00336 CALL get_qs_env(qs_env, input=input, error=error) 00337 xc_section => section_vals_get_subs_vals(input,"DFT%XC",error=error) 00338 CALL create_admm_xc_section(qs_env, xc_section, admm_env, error) 00339 admm_env%method_id = admm_control%method_id 00340 admm_env%purification_method = admm_control%purification_method 00341 admm_env%block_purification_method = admm_control%block_purification_method 00342 admm_env%block_projection_method = admm_control%block_projection_method 00343 IF( admm_env%method_id == do_admm_block_density_matrix) THEN 00344 admm_block_section => section_vals_get_subs_vals(input,& 00345 "DFT%AUXILIARY_DENSITY_MATRIX_METHOD%BLOCK_DENSITY_MATRIX_METHOD",error=error) 00346 CALL get_qs_env(qs_env=qs_env,& 00347 particle_set=particle_set,& 00348 error=error) 00349 natom = SIZE(particle_set,1) 00350 00351 CALL admm_create_block_list(admm_block_section, admm_env, natom, error) 00352 END IF 00353 END IF 00354 00355 SELECT CASE(admm_env%method_id) 00356 CASE (do_admm_basis_set_projection) 00357 SELECT CASE(admm_env%purification_method) 00358 CASE(do_admm_purify_none) 00359 CALL fit_mo_coeffs_none(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, & 00360 mos, mos_aux_fit, geometry_did_change, error) 00361 CASE(do_admm_purify_cauchy) 00362 CALL fit_mo_coeffs_none(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, & 00363 mos, mos_aux_fit, geometry_did_change, error) 00364 CASE(do_admm_purify_cauchy_subspace) 00365 CALL fit_mo_coeffs_no_diag(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, & 00366 mos, mos_aux_fit, geometry_did_change, error) 00367 CASE(do_admm_purify_mo_diag) 00368 CALL fit_mo_coeffs_diag(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, & 00369 mos, mos_aux_fit, geometry_did_change, error) 00370 CASE(do_admm_purify_mo_no_diag) 00371 CALL fit_mo_coeffs_no_diag(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, & 00372 mos, mos_aux_fit, geometry_did_change, error) 00373 END SELECT 00374 00375 CASE(do_admm_block_density_matrix) 00376 SELECT CASE(admm_env%block_purification_method) 00377 CASE(do_admm_block_purify_full) 00378 SELECT CASE(admm_env%block_projection_method) 00379 CASE (do_admm_block_aux_basis_off) 00380 CALL fit_mo_coeffs_none(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, & 00381 mos, mos_aux_fit, geometry_did_change, error) 00382 CASE (do_admm_block_aux_basis_on) 00383 END SELECT 00384 CASE(do_admm_block_purify_off) 00385 SELECT CASE(admm_env%block_projection_method) 00386 CASE (do_admm_block_aux_basis_off) 00387 CALL fit_mo_coeffs_none(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, & 00388 mos, mos_aux_fit, geometry_did_change, error) 00389 CASE (do_admm_block_aux_basis_on) 00390 END SELECT 00391 CASE(do_admm_block_purify_blocked) 00392 SELECT CASE(admm_env%block_projection_method) 00393 CASE (do_admm_block_aux_basis_off) 00394 CALL fit_mo_coeffs_blocked(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, & 00395 mos, mos_aux_fit, geometry_did_change, error) 00396 CASE (do_admm_block_aux_basis_on) 00397 END SELECT 00398 00399 END SELECT 00400 00401 END SELECT 00402 00403 CALL timestop(handle) 00404 00405 END SUBROUTINE admm_fit_mo_coeffs 00406 00407 ! ***************************************************************************** 00424 SUBROUTINE fit_mo_coeffs_diag(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, & 00425 mos, mos_aux_fit, geometry_did_change, error) 00426 00427 TYPE(qs_environment_type), POINTER :: qs_env 00428 TYPE(admm_type), POINTER :: admm_env 00429 TYPE(cp_para_env_type), POINTER :: para_env 00430 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00431 POINTER :: matrix_s_aux_fit, 00432 matrix_s_mixed 00433 TYPE(mo_set_p_type), DIMENSION(:), 00434 POINTER :: mos, mos_aux_fit 00435 LOGICAL, INTENT(IN) :: geometry_did_change 00436 TYPE(cp_error_type), INTENT(inout) :: error 00437 00438 CHARACTER(LEN=*), PARAMETER :: routineN = 'fit_mo_coeffs_diag', 00439 routineP = moduleN//':'//routineN 00440 00441 INTEGER :: handle, i, ispin, istat, 00442 nao_aux_fit, nao_orb, nmo, 00443 nspins 00444 LOGICAL :: failure 00445 REAL(dp), ALLOCATABLE, DIMENSION(:) :: eig_work 00446 TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit 00447 TYPE(section_vals_type), POINTER :: input, xc_section 00448 00449 CALL timeset(routineN,handle) 00450 00451 IF (.NOT.(ASSOCIATED(admm_env) )) THEN 00452 CALL admm_env_create(mos, mos_aux_fit, & 00453 para_env, admm_env,& 00454 error) 00455 CALL get_qs_env(qs_env, input=input, error=error) 00456 xc_section => section_vals_get_subs_vals(input,"DFT%XC",error=error) 00457 CALL create_admm_xc_section(qs_env, xc_section, admm_env, error) 00458 END IF 00459 00460 nao_aux_fit = admm_env%nao_aux_fit 00461 nao_orb = admm_env%nao_orb 00462 nspins = SIZE(mos) 00463 00464 00465 ! *** This part only depends on overlap matrices ==> needs only to be calculated if the geometry changed 00466 IF( geometry_did_change ) THEN 00467 CALL copy_dbcsr_to_fm(matrix_s_aux_fit(1)%matrix,admm_env%S_inv,error) 00468 CALL cp_fm_upper_to_full(admm_env%S_inv,admm_env%work_aux_aux,error=error) 00469 CALL cp_fm_to_fm(admm_env%S_inv, admm_env%S, error=error) 00470 00471 CALL copy_dbcsr_to_fm(matrix_s_mixed(1)%matrix,admm_env%Q,error) 00472 00473 !! Calculate S'_inverse 00474 CALL cp_fm_cholesky_decompose(admm_env%S_inv,error=error) 00475 CALL cp_fm_cholesky_invert(admm_env%S_inv,error=error) 00476 !! Symmetrize the guy 00477 CALL cp_fm_upper_to_full(admm_env%S_inv,admm_env%work_aux_aux,error=error) 00478 !! Calculate A=S'^(-1)*P 00479 CALL cp_fm_gemm('N','N',nao_aux_fit,nao_orb,nao_aux_fit,& 00480 1.0_dp,admm_env%S_inv,admm_env%Q,0.0_dp,& 00481 admm_env%A,error) 00482 00483 !! B=Q^(T)*A 00484 CALL cp_fm_gemm('T','N',nao_orb,nao_orb,nao_aux_fit,& 00485 1.0_dp,admm_env%Q,admm_env%A,0.0_dp,& 00486 admm_env%B,error) 00487 00488 END IF 00489 00490 ! *** Calculate the mo_coeffs for the fitting basis 00491 DO ispin=1,nspins 00492 nmo = admm_env%nmo(ispin) 00493 IF( nmo == 0 ) CYCLE 00494 !! Lambda = C^(T)*B*C 00495 CALL get_mo_set(mos(ispin)%mo_set,mo_coeff=mo_coeff) 00496 CALL get_mo_set(mos_aux_fit(ispin)%mo_set,mo_coeff=mo_coeff_aux_fit) 00497 CALL cp_fm_gemm('N','N',nao_orb,nmo,nao_orb,& 00498 1.0_dp,admm_env%B,mo_coeff,0.0_dp,& 00499 admm_env%work_orb_nmo(ispin)%matrix,error) 00500 CALL cp_fm_gemm('T','N',nmo,nmo,nao_orb,& 00501 1.0_dp,mo_coeff,admm_env%work_orb_nmo(ispin)%matrix,0.0_dp,& 00502 admm_env%lambda(ispin)%matrix,error) 00503 CALL cp_fm_to_fm(admm_env%lambda(ispin)%matrix, admm_env%work_nmo_nmo1(ispin)%matrix, error=error) 00504 CALL cp_fm_syevd(admm_env%work_nmo_nmo1(ispin)%matrix,admm_env%R(ispin)%matrix,& 00505 admm_env%eigvals_lambda(ispin)%eigvals%data,error=error) 00506 ALLOCATE(eig_work(nmo), STAT=istat) 00507 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00508 DO i=1,nmo 00509 eig_work(i) = 1.0_dp/SQRT(admm_env%eigvals_lambda(ispin)%eigvals%data(i)) 00510 END DO 00511 CALL cp_fm_to_fm(admm_env%R(ispin)%matrix, admm_env%work_nmo_nmo1(ispin)%matrix, error=error) 00512 CALL cp_fm_column_scale(admm_env%work_nmo_nmo1(ispin)%matrix,eig_work) 00513 CALL cp_fm_gemm('N','T',nmo,nmo,nmo,& 00514 1.0_dp,admm_env%work_nmo_nmo1(ispin)%matrix,admm_env%R(ispin)%matrix,0.0_dp,& 00515 admm_env%lambda_inv_sqrt(ispin)%matrix,error) 00516 CALL cp_fm_gemm('N','N',nao_orb,nmo,nmo,& 00517 1.0_dp,mo_coeff,admm_env%lambda_inv_sqrt(ispin)%matrix,0.0_dp,& 00518 admm_env%work_orb_nmo(ispin)%matrix,error) 00519 CALL cp_fm_gemm('N','N',nao_aux_fit,nmo,nao_orb,& 00520 1.0_dp,admm_env%A,admm_env%work_orb_nmo(ispin)%matrix, 0.0_dp,& 00521 mo_coeff_aux_fit,error) 00522 DEALLOCATE(eig_work, STAT=istat) 00523 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00524 00525 END DO 00526 00527 CALL timestop(handle) 00528 00529 END SUBROUTINE fit_mo_coeffs_diag 00530 00531 SUBROUTINE fit_mo_coeffs_no_diag(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, & 00532 mos, mos_aux_fit, geometry_did_change, error) 00533 TYPE(qs_environment_type), POINTER :: qs_env 00534 TYPE(admm_type), POINTER :: admm_env 00535 TYPE(cp_para_env_type), POINTER :: para_env 00536 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00537 POINTER :: matrix_s_aux_fit, 00538 matrix_s_mixed 00539 TYPE(mo_set_p_type), DIMENSION(:), 00540 POINTER :: mos, mos_aux_fit 00541 LOGICAL, INTENT(IN) :: geometry_did_change 00542 TYPE(cp_error_type), INTENT(inout) :: error 00543 00544 CHARACTER(LEN=*), PARAMETER :: routineN = 'fit_mo_coeffs_no_diag', 00545 routineP = moduleN//':'//routineN 00546 00547 INTEGER :: handle, ispin, nao_aux_fit, 00548 nao_orb, nmo, nspins 00549 TYPE(cp_fm_type), POINTER :: mo_coeff 00550 TYPE(section_vals_type), POINTER :: input, xc_section 00551 00552 CALL timeset(routineN,handle) 00553 00554 IF (.NOT.(ASSOCIATED(admm_env) )) THEN 00555 CALL admm_env_create(mos, mos_aux_fit, & 00556 para_env, admm_env,& 00557 error) 00558 CALL get_qs_env(qs_env, input=input, error=error) 00559 xc_section => section_vals_get_subs_vals(input,"DFT%XC",error=error) 00560 CALL create_admm_xc_section(qs_env, xc_section, admm_env, error) 00561 END IF 00562 00563 nao_aux_fit = admm_env%nao_aux_fit 00564 nao_orb = admm_env%nao_orb 00565 nspins = SIZE(mos) 00566 00567 00568 ! *** This part only depends on overlap matrices ==> needs only to be calculated if the geometry changed 00569 00570 IF( geometry_did_change ) THEN 00571 CALL copy_dbcsr_to_fm(matrix_s_aux_fit(1)%matrix,admm_env%S_inv,error) 00572 CALL cp_fm_upper_to_full(admm_env%S_inv,admm_env%work_aux_aux,error=error) 00573 CALL cp_fm_to_fm(admm_env%S_inv, admm_env%S, error=error) 00574 00575 CALL copy_dbcsr_to_fm(matrix_s_mixed(1)%matrix,admm_env%Q,error) 00576 00577 !! Calculate S'_inverse 00578 CALL cp_fm_cholesky_decompose(admm_env%S_inv,error=error) 00579 CALL cp_fm_cholesky_invert(admm_env%S_inv,error=error) 00580 !! Symmetrize the guy 00581 CALL cp_fm_upper_to_full(admm_env%S_inv,admm_env%work_aux_aux,error=error) 00582 !! Calculate A=S'^(-1)*Q 00583 CALL cp_fm_gemm('N','N',nao_aux_fit,nao_orb,nao_aux_fit,& 00584 1.0_dp,admm_env%S_inv,admm_env%Q,0.0_dp,& 00585 admm_env%A,error) 00586 00587 !! B=Q^(T)*A 00588 CALL cp_fm_gemm('T','N',nao_orb,nao_orb,nao_aux_fit,& 00589 1.0_dp,admm_env%Q,admm_env%A,0.0_dp,& 00590 admm_env%B,error) 00591 END IF 00592 00593 DO ispin = 1,nspins 00594 nmo = admm_env%nmo(ispin) 00595 CALL get_mo_set(mos(ispin)%mo_set,mo_coeff=mo_coeff) 00596 00597 ! Calculate Lambda inverse for later usage 00598 CALL cp_fm_gemm('N','N',nao_orb,nmo,nao_orb,& 00599 1.0_dp,admm_env%B,mo_coeff,0.0_dp,& 00600 admm_env%work_orb_nmo(ispin)%matrix,error) 00601 CALL cp_fm_gemm('T','N',nmo,nmo,nao_orb,& 00602 1.0_dp,mo_coeff,admm_env%work_orb_nmo(ispin)%matrix,0.0_dp,& 00603 admm_env%lambda(ispin)%matrix,error) 00604 CALL cp_fm_to_fm(admm_env%lambda(ispin)%matrix, admm_env%work_nmo_nmo1(ispin)%matrix, error=error) 00605 00606 CALL cp_fm_cholesky_decompose(admm_env%work_nmo_nmo1(ispin)%matrix,error=error) 00607 CALL cp_fm_cholesky_invert(admm_env%work_nmo_nmo1(ispin)%matrix,error=error) 00608 !! Symmetrize the guy 00609 CALL cp_fm_upper_to_full(admm_env%work_nmo_nmo1(ispin)%matrix,admm_env%lambda_inv(ispin)%matrix,error=error) 00610 CALL cp_fm_to_fm(admm_env%work_nmo_nmo1(ispin)%matrix,admm_env%lambda_inv(ispin)%matrix,error=error) 00611 00612 00613 !! ** C_hat = AC 00614 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nmo, nao_orb,& 00615 1.0_dp,admm_env%A,mo_coeff,0.0_dp,& 00616 admm_env%C_hat(ispin)%matrix,error) 00617 END DO 00618 00619 CALL timestop(handle) 00620 00621 END SUBROUTINE fit_mo_coeffs_no_diag 00622 00623 SUBROUTINE fit_mo_coeffs_none(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, & 00624 mos, mos_aux_fit, geometry_did_change, error) 00625 TYPE(qs_environment_type), POINTER :: qs_env 00626 TYPE(admm_type), POINTER :: admm_env 00627 TYPE(cp_para_env_type), POINTER :: para_env 00628 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00629 POINTER :: matrix_s_aux_fit, 00630 matrix_s_mixed 00631 TYPE(mo_set_p_type), DIMENSION(:), 00632 POINTER :: mos, mos_aux_fit 00633 LOGICAL, INTENT(IN) :: geometry_did_change 00634 TYPE(cp_error_type), INTENT(inout) :: error 00635 00636 CHARACTER(LEN=*), PARAMETER :: routineN = 'fit_mo_coeffs_none', 00637 routineP = moduleN//':'//routineN 00638 00639 INTEGER :: handle, ispin, nao_aux_fit, 00640 nao_orb, nmo, nmo_mos, nspins 00641 REAL(KIND=dp), DIMENSION(:), POINTER :: occ_num, occ_num_aux 00642 TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit 00643 TYPE(section_vals_type), POINTER :: input, xc_section 00644 00645 CALL timeset(routineN,handle) 00646 00647 IF (.NOT.(ASSOCIATED(admm_env) )) THEN 00648 CALL admm_env_create(mos, mos_aux_fit, & 00649 para_env, admm_env,& 00650 error) 00651 CALL get_qs_env(qs_env, input=input, error=error) 00652 xc_section => section_vals_get_subs_vals(input,"DFT%XC",error=error) 00653 CALL create_admm_xc_section(qs_env, xc_section, admm_env, error) 00654 END IF 00655 00656 nao_aux_fit = admm_env%nao_aux_fit 00657 nao_orb = admm_env%nao_orb 00658 nspins = SIZE(mos) 00659 00660 00661 ! *** This part only depends on overlap matrices ==> needs only to be calculated if the geometry changed 00662 00663 IF( geometry_did_change ) THEN 00664 CALL copy_dbcsr_to_fm(matrix_s_aux_fit(1)%matrix,admm_env%S_inv,error) 00665 CALL cp_fm_upper_to_full(admm_env%S_inv,admm_env%work_aux_aux,error=error) 00666 CALL cp_fm_to_fm(admm_env%S_inv, admm_env%S, error=error) 00667 00668 CALL copy_dbcsr_to_fm(matrix_s_mixed(1)%matrix,admm_env%Q,error) 00669 00670 !! Calculate S'_inverse 00671 CALL cp_fm_cholesky_decompose(admm_env%S_inv,error=error) 00672 CALL cp_fm_cholesky_invert(admm_env%S_inv,error=error) 00673 !! Symmetrize the guy 00674 CALL cp_fm_upper_to_full(admm_env%S_inv,admm_env%work_aux_aux,error=error) 00675 !! Calculate A=S'^(-1)*Q 00676 CALL cp_fm_gemm('N','N',nao_aux_fit,nao_orb,nao_aux_fit,& 00677 1.0_dp,admm_env%S_inv,admm_env%Q,0.0_dp,& 00678 admm_env%A,error) 00679 END IF 00680 00681 DO ispin = 1,nspins 00682 nmo = admm_env%nmo(ispin) 00683 CALL get_mo_set(mos(ispin)%mo_set,mo_coeff=mo_coeff, occupation_numbers=occ_num,nmo=nmo_mos) 00684 CALL get_mo_set(mos_aux_fit(ispin)%mo_set,mo_coeff=mo_coeff_aux_fit,& 00685 occupation_numbers=occ_num_aux) 00686 00687 CALL cp_fm_gemm('N','N',nao_aux_fit,nmo,nao_orb,& 00688 1.0_dp,admm_env%A,mo_coeff,0.0_dp,& 00689 mo_coeff_aux_fit,error) 00690 00691 occ_num_aux(1:nmo) = occ_num(1:nmo) 00692 ! XXXX should only be done first time XXXX 00693 CALL cp_fm_set_all(admm_env%lambda(ispin)%matrix,0.0_dp,1.0_dp,error) 00694 CALL cp_fm_set_all(admm_env%lambda_inv(ispin)%matrix,0.0_dp,1.0_dp,error) 00695 CALL cp_fm_set_all(admm_env%lambda_inv_sqrt(ispin)%matrix,0.0_dp,1.0_dp,error) 00696 END DO 00697 00698 CALL timestop(handle) 00699 00700 END SUBROUTINE fit_mo_coeffs_none 00701 00702 SUBROUTINE fit_mo_coeffs_blocked(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, & 00703 mos, mos_aux_fit, geometry_did_change, error) 00704 TYPE(qs_environment_type), POINTER :: qs_env 00705 TYPE(admm_type), POINTER :: admm_env 00706 TYPE(cp_para_env_type), POINTER :: para_env 00707 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00708 POINTER :: matrix_s_aux_fit, 00709 matrix_s_mixed 00710 TYPE(mo_set_p_type), DIMENSION(:), 00711 POINTER :: mos, mos_aux_fit 00712 LOGICAL, INTENT(IN) :: geometry_did_change 00713 TYPE(cp_error_type), INTENT(inout) :: error 00714 00715 CHARACTER(LEN=*), PARAMETER :: routineN = 'fit_mo_coeffs_blocked', 00716 routineP = moduleN//':'//routineN 00717 00718 INTEGER :: blk, handle, iatom, jatom, 00719 nao_aux_fit, nao_orb, nspins 00720 REAL(dp), DIMENSION(:, :), POINTER :: sparse_block 00721 TYPE(cp_dbcsr_iterator) :: iter 00722 TYPE(cp_dbcsr_type), POINTER :: matrix_s_tilde 00723 TYPE(section_vals_type), POINTER :: input, xc_section 00724 00725 CALL timeset(routineN,handle) 00726 00727 IF (.NOT.(ASSOCIATED(admm_env) )) THEN 00728 CALL admm_env_create(mos, mos_aux_fit, & 00729 para_env, admm_env,& 00730 error) 00731 CALL get_qs_env(qs_env, input=input, error=error) 00732 xc_section => section_vals_get_subs_vals(input,"DFT%XC",error=error) 00733 CALL create_admm_xc_section(qs_env, xc_section, admm_env, error) 00734 END IF 00735 00736 nao_aux_fit = admm_env%nao_aux_fit 00737 nao_orb = admm_env%nao_orb 00738 nspins = SIZE(mos) 00739 00740 00741 ! *** This part only depends on overlap matrices ==> needs only to be calculated if the geometry changed 00742 00743 IF( geometry_did_change ) THEN 00744 NULLIFY(matrix_s_tilde) 00745 ALLOCATE(matrix_s_tilde) 00746 CALL cp_dbcsr_init (matrix_s_tilde, error) 00747 CALL cp_dbcsr_create(matrix_s_tilde, 'MATRIX s_tilde', & 00748 cp_dbcsr_distribution(matrix_s_aux_fit(1)%matrix), dbcsr_type_symmetric, & 00749 cp_dbcsr_row_block_sizes(matrix_s_aux_fit(1)%matrix),& 00750 cp_dbcsr_col_block_sizes(matrix_s_aux_fit(1)%matrix), & 00751 cp_dbcsr_get_num_blocks(matrix_s_aux_fit(1)%matrix), & 00752 cp_dbcsr_get_data_size(matrix_s_aux_fit(1)%matrix),& 00753 cp_dbcsr_get_data_type(matrix_s_aux_fit(1)%matrix), & 00754 error=error) 00755 00756 CALL cp_dbcsr_copy(matrix_s_tilde, matrix_s_aux_fit(1)%matrix, error=error) 00757 00758 CALL cp_dbcsr_iterator_start(iter, matrix_s_tilde) 00759 DO WHILE (cp_dbcsr_iterator_blocks_left (iter)) 00760 CALL cp_dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block, blk) 00761 IF( admm_env%block_map(iatom,jatom) == 0 ) THEN 00762 sparse_block = 0.0_dp 00763 END IF 00764 END DO 00765 CALL cp_dbcsr_iterator_stop(iter) 00766 00767 00768 CALL copy_dbcsr_to_fm(matrix_s_tilde,admm_env%S_inv,error) 00769 00770 CALL cp_dbcsr_deallocate_matrix(matrix_s_tilde,error) 00771 00772 CALL cp_fm_upper_to_full(admm_env%S_inv,admm_env%work_aux_aux,error=error) 00773 CALL cp_fm_to_fm(admm_env%S_inv, admm_env%S, error=error) 00774 00775 CALL copy_dbcsr_to_fm(matrix_s_mixed(1)%matrix,admm_env%Q,error) 00776 00777 !! Calculate S'_inverse 00778 CALL cp_fm_cholesky_decompose(admm_env%S_inv,error=error) 00779 CALL cp_fm_cholesky_invert(admm_env%S_inv,error=error) 00780 !! Symmetrize the guy 00781 CALL cp_fm_upper_to_full(admm_env%S_inv,admm_env%work_aux_aux,error=error) 00782 !! Calculate A=S'^(-1)*Q 00783 CALL cp_fm_set_all(admm_env%A, 0.0_dp, 1.0_dp, error) 00784 END IF 00785 00786 CALL timestop(handle) 00787 00788 END SUBROUTINE fit_mo_coeffs_blocked 00789 00790 00791 ! ***************************************************************************** 00809 SUBROUTINE merge_mo_derivs_diag(ispin, admm_env, mo_set, mo_coeff, mo_coeff_aux_fit, mo_derivs, & 00810 mo_derivs_aux_fit, matrix_ks_aux_fit, error) 00811 INTEGER, INTENT(IN) :: ispin 00812 TYPE(admm_type), POINTER :: admm_env 00813 TYPE(mo_set_type), POINTER :: mo_set 00814 TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit 00815 TYPE(cp_fm_p_type), DIMENSION(:), 00816 POINTER :: mo_derivs, mo_derivs_aux_fit 00817 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00818 POINTER :: matrix_ks_aux_fit 00819 TYPE(cp_error_type), INTENT(INOUT) :: error 00820 00821 CHARACTER(LEN=*), PARAMETER :: routineN = 'merge_mo_derivs_diag', 00822 routineP = moduleN//':'//routineN 00823 00824 INTEGER :: handle, i, j, nao_aux_fit, 00825 nao_orb, nmo, stat 00826 INTEGER, SAVE :: counter = 0 00827 LOGICAL :: failure 00828 REAL(dp) :: eig_diff, pole, tmp32, tmp52, 00829 tmp72, tmp92 00830 REAL(dp), DIMENSION(:), POINTER :: occupation_numbers, 00831 scaling_factor 00832 00833 failure = .FALSE. 00834 00835 CALL timeset(routineN,handle) 00836 00837 counter = counter + 1 00838 00839 nao_aux_fit = admm_env%nao_aux_fit 00840 nao_orb = admm_env%nao_orb 00841 nmo = admm_env%nmo(ispin) 00842 00843 CALL copy_dbcsr_to_fm(matrix_ks_aux_fit(ispin)%matrix,admm_env%K(ispin)%matrix,error) 00844 CALL cp_fm_upper_to_full(admm_env%K(ispin)%matrix,admm_env%work_aux_aux,error=error) 00845 00846 CALL cp_fm_gemm('N','N', nao_aux_fit, nmo, nao_aux_fit,& 00847 1.0_dp,admm_env%K(ispin)%matrix,mo_coeff_aux_fit,0.0_dp,& 00848 admm_env%H(ispin)%matrix,error) 00849 00850 CALL get_mo_set(mo_set=mo_set, occupation_numbers=occupation_numbers) 00851 ALLOCATE(scaling_factor(SIZE(occupation_numbers)),stat=stat) 00852 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00853 scaling_factor = 2.0_dp*occupation_numbers 00854 00855 CALL cp_fm_column_scale(admm_env%H(ispin)%matrix,scaling_factor) 00856 00857 CALL cp_fm_to_fm(admm_env%H(ispin)%matrix, mo_derivs_aux_fit(ispin)%matrix, error=error) 00858 00859 ! *** Add first term 00860 CALL cp_fm_gemm('N','T', nao_aux_fit, nmo, nmo,& 00861 1.0_dp,admm_env%H(ispin)%matrix,admm_env%lambda_inv_sqrt(ispin)%matrix,0.0_dp,& 00862 admm_env%work_aux_nmo(ispin)%matrix,error) 00863 CALL cp_fm_gemm('T','N', nao_orb, nmo, nao_aux_fit,& 00864 1.0_dp,admm_env%A,admm_env%work_aux_nmo(ispin)%matrix,0.0_dp,& 00865 admm_env%mo_derivs_tmp(ispin)%matrix,error) 00866 00867 00868 ! *** Construct Matrix M for Hadamard Product 00869 pole = 0.0_dp 00870 DO i=1,nmo 00871 DO j=i,nmo 00872 eig_diff = ( admm_env%eigvals_lambda(ispin)%eigvals%data(i) -& 00873 admm_env%eigvals_lambda(ispin)%eigvals%data(j) ) 00874 ! *** two eigenvalues could be the degenerated. In that case use 2nd order formula for the poles 00875 IF( ABS(eig_diff) < 0.0001_dp ) THEN 00876 tmp32 = 1.0_dp/SQRT(admm_env%eigvals_lambda(ispin)%eigvals%data(j))**3 00877 tmp52 = tmp32/admm_env%eigvals_lambda(ispin)%eigvals%data(j)*eig_diff 00878 tmp72 = tmp52/admm_env%eigvals_lambda(ispin)%eigvals%data(j)*eig_diff 00879 tmp92 = tmp72/admm_env%eigvals_lambda(ispin)%eigvals%data(j)*eig_diff 00880 00881 pole = -0.5_dp*tmp32 + 3.0_dp/8.0_dp*tmp52 - 5.0_dp/16.0_dp*tmp72 + 35.0_dp/128.0_dp*tmp92 00882 CALL cp_fm_set_element(admm_env%M(ispin)%matrix,i,j,pole,error) 00883 ELSE 00884 pole = 1.0_dp/SQRT(admm_env%eigvals_lambda(ispin)%eigvals%data(i)) 00885 pole = pole - 1.0_dp/SQRT(admm_env%eigvals_lambda(ispin)%eigvals%data(j)) 00886 pole = pole/(admm_env%eigvals_lambda(ispin)%eigvals%data(i)-& 00887 admm_env%eigvals_lambda(ispin)%eigvals%data(j)) 00888 CALL cp_fm_set_element(admm_env%M(ispin)%matrix,i,j,pole,error) 00889 END IF 00890 END DO 00891 END DO 00892 CALL cp_fm_upper_to_full(admm_env%M(ispin)%matrix,admm_env%work_nmo_nmo1(ispin)%matrix,error=error) 00893 00894 ! *** 2nd term to be added to fm_H 00895 00896 !! Part 1: B^(T)*C* R*[R^(T)*c^(T)*A^(T)*H_aux_fit*R x M]*R^(T) 00897 !! Part 2: B*C*(R*[R^(T)*c^(T)*A^(T)*H_aux_fit*R x M]*R^(T))^(T) 00898 00899 ! *** H'*R 00900 CALL cp_fm_gemm('N','N', nao_aux_fit, nmo, nmo,& 00901 1.0_dp,admm_env%H(ispin)%matrix,admm_env%R(ispin)%matrix,0.0_dp,& 00902 admm_env%work_aux_nmo(ispin)%matrix,error) 00903 ! *** A^(T)*H'*R 00904 CALL cp_fm_gemm('T','N', nao_orb, nmo, nao_aux_fit,& 00905 1.0_dp,admm_env%A,admm_env%work_aux_nmo(ispin)%matrix,0.0_dp,& 00906 admm_env%work_orb_nmo(ispin)%matrix,error) 00907 ! *** c^(T)*A^(T)*H'*R 00908 CALL cp_fm_gemm('T','N', nmo, nmo, nao_orb,& 00909 1.0_dp,mo_coeff,admm_env%work_orb_nmo(ispin)%matrix,0.0_dp,& 00910 admm_env%work_nmo_nmo1(ispin)%matrix,error) 00911 ! *** R^(T)*c^(T)*A^(T)*H'*R 00912 CALL cp_fm_gemm('T','N', nmo, nmo, nmo,& 00913 1.0_dp,admm_env%R(ispin)%matrix,admm_env%work_nmo_nmo1(ispin)%matrix,0.0_dp,& 00914 admm_env%work_nmo_nmo2(ispin)%matrix,error) 00915 ! *** R^(T)*c^(T)*A^(T)*H'*R x M 00916 CALL cp_fm_schur_product(admm_env%work_nmo_nmo2(ispin)%matrix,& 00917 admm_env%M(ispin)%matrix,admm_env%work_nmo_nmo1(ispin)%matrix,error) 00918 ! *** R* (R^(T)*c^(T)*A^(T)*H'*R x M) 00919 CALL cp_fm_gemm('N','N', nmo, nmo, nmo,& 00920 1.0_dp,admm_env%R(ispin)%matrix,admm_env%work_nmo_nmo1(ispin)%matrix,0.0_dp,& 00921 admm_env%work_nmo_nmo2(ispin)%matrix,error) 00922 00923 ! *** R* (R^(T)*c^(T)*A^(T)*H'*R x M) *R^(T) 00924 CALL cp_fm_gemm('N','T', nmo, nmo, nmo,& 00925 1.0_dp,admm_env%work_nmo_nmo2(ispin)%matrix,admm_env%R(ispin)%matrix,0.0_dp,& 00926 admm_env%R_schur_R_t(ispin)%matrix,error) 00927 00928 ! *** B^(T)*c 00929 CALL cp_fm_gemm('T','N', nao_orb, nmo, nao_orb,& 00930 1.0_dp,admm_env%B,mo_coeff,0.0_dp,& 00931 admm_env%work_orb_nmo(ispin)%matrix,error) 00932 00933 ! *** Add first term to fm_H 00934 ! *** B^(T)*c* R* (R^(T)*c^(T)*A^(T)*H'*R x M) *R^(T) 00935 CALL cp_fm_gemm('N','N', nao_orb, nmo, nmo,& 00936 1.0_dp,admm_env%work_orb_nmo(ispin)%matrix,admm_env%R_schur_R_t(ispin)%matrix,1.0_dp,& 00937 admm_env%mo_derivs_tmp(ispin)%matrix,error) 00938 00939 ! *** Add second term to fm_H 00940 ! *** B*C *[ R* (R^(T)*c^(T)*A^(T)*H'*R x M) *R^(T)]^(T) 00941 CALL cp_fm_gemm('N','T', nao_orb, nmo, nmo,& 00942 1.0_dp,admm_env%work_orb_nmo(ispin)%matrix,admm_env%R_schur_R_t(ispin)%matrix,1.0_dp,& 00943 admm_env%mo_derivs_tmp(ispin)%matrix,error) 00944 00945 DO i = 1,SIZE(scaling_factor) 00946 scaling_factor(i) = 1.0_dp/scaling_factor(i) 00947 END DO 00948 00949 CALL cp_fm_column_scale(admm_env%mo_derivs_tmp(ispin)%matrix,scaling_factor) 00950 00951 CALL cp_fm_scale_and_add(1.0_dp,mo_derivs(ispin)%matrix,1.0_dp,admm_env%mo_derivs_tmp(ispin)%matrix,error) 00952 00953 DEALLOCATE(scaling_factor, stat=stat) 00954 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00955 00956 CALL timestop(handle) 00957 00958 END SUBROUTINE merge_mo_derivs_diag 00959 00960 SUBROUTINE calc_aux_mo_derivs_none(ispin, admm_env, mo_set, mo_coeff, mo_coeff_aux_fit, & 00961 matrix_ks_aux_fit, error) 00962 INTEGER, INTENT(IN) :: ispin 00963 TYPE(admm_type), POINTER :: admm_env 00964 TYPE(mo_set_type), POINTER :: mo_set 00965 TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit 00966 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00967 POINTER :: matrix_ks_aux_fit 00968 TYPE(cp_error_type), INTENT(INOUT) :: error 00969 00970 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_aux_mo_derivs_none', 00971 routineP = moduleN//':'//routineN 00972 00973 INTEGER :: handle, nao_aux_fit, nao_orb, 00974 nmo, stat 00975 INTEGER, SAVE :: counter = 0 00976 LOGICAL :: failure 00977 REAL(dp), DIMENSION(:), POINTER :: occupation_numbers, 00978 scaling_factor 00979 00980 failure = .FALSE. 00981 00982 CALL timeset(routineN,handle) 00983 00984 counter = counter + 1 00985 00986 nao_aux_fit = admm_env%nao_aux_fit 00987 nao_orb = admm_env%nao_orb 00988 nmo = admm_env%nmo(ispin) 00989 00990 ! just calculate the mo derivs in the aux basis 00991 ! only needs to be done on the converged ks matrix for the force calc 00992 ! Note with OT and purification NONE, the merging of the derivs 00993 ! happens implicitly because the KS matrices have been already been merged 00994 ! and adding them here would be double counting. 00995 00996 CALL copy_dbcsr_to_fm(matrix_ks_aux_fit(ispin)%matrix,admm_env%K(ispin)%matrix,error) 00997 CALL cp_fm_upper_to_full(admm_env%K(ispin)%matrix,admm_env%work_aux_aux,error=error) 00998 00999 CALL cp_fm_gemm('N','N', nao_aux_fit, nmo, nao_aux_fit,& 01000 1.0_dp,admm_env%K(ispin)%matrix,mo_coeff_aux_fit,0.0_dp,& 01001 admm_env%H(ispin)%matrix,error) 01002 01003 CALL get_mo_set(mo_set=mo_set, occupation_numbers=occupation_numbers) 01004 ALLOCATE(scaling_factor(SIZE(occupation_numbers)),stat=stat) 01005 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01006 scaling_factor = 2.0_dp*occupation_numbers 01007 CALL cp_fm_column_scale(admm_env%H(ispin)%matrix,scaling_factor) 01008 01009 DEALLOCATE(scaling_factor, stat=stat) 01010 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01011 01012 CALL timestop(handle) 01013 01014 END SUBROUTINE calc_aux_mo_derivs_none 01015 01016 SUBROUTINE merge_mo_derivs_no_diag(ispin, admm_env, mo_set, mo_coeff, mo_coeff_aux_fit, mo_derivs, & 01017 mo_derivs_aux_fit, matrix_ks_aux_fit, error) 01018 INTEGER, INTENT(IN) :: ispin 01019 TYPE(admm_type), POINTER :: admm_env 01020 TYPE(mo_set_type), POINTER :: mo_set 01021 TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit 01022 TYPE(cp_fm_p_type), DIMENSION(:), 01023 POINTER :: mo_derivs, mo_derivs_aux_fit 01024 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01025 POINTER :: matrix_ks_aux_fit 01026 TYPE(cp_error_type), INTENT(INOUT) :: error 01027 01028 CHARACTER(LEN=*), PARAMETER :: routineN = 'merge_mo_derivs_no_diag', 01029 routineP = moduleN//':'//routineN 01030 01031 INTEGER :: handle, nao_aux_fit, nao_orb, 01032 nmo, stat 01033 INTEGER, SAVE :: counter = 0 01034 LOGICAL :: failure 01035 REAL(dp), DIMENSION(:), POINTER :: occupation_numbers, 01036 scaling_factor 01037 01038 failure = .FALSE. 01039 01040 CALL timeset(routineN,handle) 01041 01042 counter = counter + 1 01043 01044 nao_aux_fit = admm_env%nao_aux_fit 01045 nao_orb = admm_env%nao_orb 01046 nmo = admm_env%nmo(ispin) 01047 01048 CALL copy_dbcsr_to_fm(matrix_ks_aux_fit(ispin)%matrix,admm_env%K(ispin)%matrix,error) 01049 CALL cp_fm_upper_to_full(admm_env%K(ispin)%matrix,admm_env%work_aux_aux,error=error) 01050 01051 CALL get_mo_set(mo_set=mo_set, occupation_numbers=occupation_numbers) 01052 ALLOCATE(scaling_factor(SIZE(occupation_numbers)),stat=stat) 01053 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01054 scaling_factor = 0.5_dp 01055 01056 01057 !! ** calculate first part 01058 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nmo, nmo,& 01059 1.0_dp,admm_env%C_hat(ispin)%matrix,admm_env%lambda_inv(ispin)%matrix,0.0_dp,& 01060 admm_env%work_aux_nmo(ispin)%matrix,error) 01061 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nmo, nao_aux_fit,& 01062 1.0_dp,admm_env%K(ispin)%matrix,admm_env%work_aux_nmo(ispin)%matrix,0.0_dp,& 01063 admm_env%work_aux_nmo2(ispin)%matrix,error) 01064 CALL cp_fm_gemm('T', 'N', nao_orb, nmo, nao_aux_fit,& 01065 2.0_dp,admm_env%A,admm_env%work_aux_nmo2(ispin)%matrix,0.0_dp,& 01066 admm_env%mo_derivs_tmp(ispin)%matrix,error) 01067 !! ** calculate second part 01068 CALL cp_fm_gemm('T', 'N', nmo, nmo, nao_aux_fit,& 01069 1.0_dp,admm_env%work_aux_nmo(ispin)%matrix,admm_env%work_aux_nmo2(ispin)%matrix,0.0_dp,& 01070 admm_env%work_orb_orb,error) 01071 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nmo, nmo,& 01072 1.0_dp,admm_env%C_hat(ispin)%matrix,admm_env%work_orb_orb,0.0_dp,& 01073 admm_env%work_aux_orb,error) 01074 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nmo, nao_aux_fit,& 01075 1.0_dp,admm_env%S,admm_env%work_aux_orb,0.0_dp,& 01076 admm_env%work_aux_nmo(ispin)%matrix,error) 01077 CALL cp_fm_gemm('T', 'N', nao_orb, nmo, nao_aux_fit,& 01078 -2.0_dp,admm_env%A,admm_env%work_aux_nmo(ispin)%matrix,1.0_dp,& 01079 admm_env%mo_derivs_tmp(ispin)%matrix,error) 01080 01081 CALL cp_fm_column_scale(admm_env%mo_derivs_tmp(ispin)%matrix,scaling_factor) 01082 01083 CALL cp_fm_scale_and_add(1.0_dp,mo_derivs(ispin)%matrix,1.0_dp,admm_env%mo_derivs_tmp(ispin)%matrix,error) 01084 01085 DEALLOCATE(scaling_factor, stat=stat) 01086 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01087 01088 CALL timestop(handle) 01089 01090 END SUBROUTINE merge_mo_derivs_no_diag 01091 01092 01093 ! ***************************************************************************** 01127 SUBROUTINE calc_mixed_overlap_force(qs_env, para_env, ispin, admm_env, mo_coeff, & 01128 matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, logger,& 01129 iw, use_virial, error) 01130 01131 TYPE(qs_environment_type), POINTER :: qs_env 01132 TYPE(cp_para_env_type), POINTER :: para_env 01133 INTEGER, INTENT(IN) :: ispin 01134 TYPE(admm_type), POINTER :: admm_env 01135 TYPE(cp_fm_type), POINTER :: mo_coeff 01136 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01137 POINTER :: matrix_s_aux_fit, 01138 matrix_s_aux_fit_vs_orb 01139 TYPE(cp_logger_type), POINTER :: logger 01140 INTEGER :: iw 01141 LOGICAL, INTENT(IN) :: use_virial 01142 TYPE(cp_error_type), INTENT(INOUT) :: error 01143 01144 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_mixed_overlap_force', 01145 routineP = moduleN//':'//routineN 01146 01147 INTEGER :: handle, nao_aux_fit, nao_orb, 01148 neighbor_list_id, nmo 01149 LOGICAL :: failure 01150 TYPE(cp_dbcsr_type), POINTER :: matrix_w_q, matrix_w_s 01151 01152 CALL timeset(routineN,handle) 01153 failure=.FALSE. 01154 01155 NULLIFY(matrix_w_q, matrix_w_s) 01156 01157 nao_aux_fit = admm_env%nao_aux_fit 01158 nao_orb = admm_env%nao_orb 01159 nmo = admm_env%nmo(ispin) 01160 01161 ! *** forces are only implemented for mo_diag or none and basis_projection *** 01162 IF (.NOT.(admm_env%method_id==do_admm_basis_set_projection .AND. & 01163 (admm_env%purification_method==do_admm_purify_mo_diag .OR. & 01164 admm_env%purification_method==do_admm_purify_none))) THEN 01165 CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) 01166 ENDIF 01167 01168 ! *** Create sparse work matrices 01169 CALL get_qs_env(qs_env=qs_env,& 01170 neighbor_list_id=neighbor_list_id,& 01171 error=error) 01172 01173 ALLOCATE(matrix_w_s) 01174 CALL cp_dbcsr_init (matrix_w_s, error) 01175 CALL cp_dbcsr_create(matrix_w_s, 'W MATRIX AUX S', & 01176 cp_dbcsr_distribution(matrix_s_aux_fit(1)%matrix), dbcsr_type_no_symmetry, & 01177 cp_dbcsr_row_block_sizes(matrix_s_aux_fit(1)%matrix),& 01178 cp_dbcsr_col_block_sizes(matrix_s_aux_fit(1)%matrix), & 01179 cp_dbcsr_get_num_blocks(matrix_s_aux_fit(1)%matrix), & 01180 cp_dbcsr_get_data_size(matrix_s_aux_fit(1)%matrix),& 01181 cp_dbcsr_get_data_type(matrix_s_aux_fit(1)%matrix), & 01182 error=error) 01183 CALL cp_dbcsr_alloc_block_from_nbl(matrix_w_s,qs_env%sab_aux_fit_asymm,error=error) 01184 01185 ALLOCATE(matrix_w_q) 01186 CALL cp_dbcsr_init(matrix_w_q, error=error) 01187 CALL cp_dbcsr_copy(matrix_w_q,matrix_s_aux_fit_vs_orb(1)%matrix,& 01188 "W MATRIX AUX Q",error=error) 01189 01190 ! *** S'^(-T)*H' 01191 IF (.NOT. admm_env%purification_method==do_admm_purify_none) THEN 01192 CALL cp_fm_gemm('T', 'N', nao_aux_fit, nmo, nao_aux_fit,& 01193 1.0_dp,admm_env%S_inv,qs_env%mo_derivs_aux_fit(ispin)%matrix,0.0_dp,& 01194 admm_env%work_aux_nmo(ispin)%matrix,error) 01195 ELSE 01196 01197 CALL cp_fm_gemm('T', 'N', nao_aux_fit, nmo, nao_aux_fit,& 01198 1.0_dp,admm_env%S_inv,admm_env%H(ispin)%matrix,0.0_dp,& 01199 admm_env%work_aux_nmo(ispin)%matrix,error) 01200 END IF 01201 01202 ! *** S'^(-T)*H'*Lambda^(-T/2) 01203 CALL cp_fm_gemm('N', 'T', nao_aux_fit,nmo, nmo,& 01204 1.0_dp,admm_env%work_aux_nmo(ispin)%matrix,admm_env%lambda_inv_sqrt(ispin)%matrix,0.0_dp,& 01205 admm_env%work_aux_nmo2(ispin)%matrix,error) 01206 01207 ! *** C*Lambda^(-1/2)*H'^(T)*S'^(-1) minus sign due to force = -dE/dR 01208 CALL cp_fm_gemm('N', 'T', nao_aux_fit, nao_orb, nmo,& 01209 -1.0_dp,admm_env%work_aux_nmo2(ispin)%matrix,mo_coeff,0.0_dp,& 01210 admm_env%work_aux_orb,error) 01211 01212 ! *** A*C*Lambda^(-1/2)*H'^(T)*S'^(-1), minus sign to recover from above 01213 CALL cp_fm_gemm('N', 'T', nao_aux_fit, nao_aux_fit, nao_orb,& 01214 -1.0_dp,admm_env%work_aux_orb,admm_env%A,0.0_dp,& 01215 admm_env%work_aux_aux,error) 01216 01217 01218 IF (.NOT. (admm_env%purification_method==do_admm_purify_none)) THEN 01219 ! *** C*Y 01220 CALL cp_fm_gemm('N', 'N', nao_orb, nmo, nmo,& 01221 1.0_dp,mo_coeff,admm_env%R_schur_R_t(ispin)%matrix,0.0_dp,& 01222 admm_env%work_orb_nmo(ispin)%matrix,error) 01223 ! *** C*Y^(T)*C^(T) 01224 CALL cp_fm_gemm('N', 'T', nao_orb, nao_orb, nmo,& 01225 1.0_dp,mo_coeff,admm_env%work_orb_nmo(ispin)%matrix,0.0_dp,& 01226 admm_env%work_orb_orb,error) 01227 ! *** A*C*Y^(T)*C^(T) Add to work aux_orb, minus sign due to force = -dE/dR 01228 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nao_orb, nao_orb,& 01229 -1.0_dp,admm_env%A,admm_env%work_orb_orb,1.0_dp,& 01230 admm_env%work_aux_orb,error) 01231 01232 ! *** C*Y^(T) 01233 CALL cp_fm_gemm('N', 'T', nao_orb, nmo, nmo,& 01234 1.0_dp,mo_coeff,admm_env%R_schur_R_t(ispin)%matrix,0.0_dp,& 01235 admm_env%work_orb_nmo(ispin)%matrix,error) 01236 ! *** C*Y*C^(T) 01237 CALL cp_fm_gemm('N', 'T', nao_orb, nao_orb, nmo,& 01238 1.0_dp,mo_coeff,admm_env%work_orb_nmo(ispin)%matrix,0.0_dp,& 01239 admm_env%work_orb_orb,error) 01240 ! *** A*C*Y*C^(T) Add to work aux_orb, minus sign due to -dE/dR 01241 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nao_orb, nao_orb,& 01242 -1.0_dp,admm_env%A,admm_env%work_orb_orb,1.0_dp,& 01243 admm_env%work_aux_orb,error) 01244 END IF 01245 01246 ! *** copy to sparse matrix 01247 CALL copy_fm_to_dbcsr(admm_env%work_aux_orb, matrix_w_q,keep_sparsity=.TRUE.,& 01248 error=error) 01249 01250 IF (.NOT. (admm_env%purification_method==do_admm_purify_none)) THEN 01251 ! *** A*C*Y^(T)*C^(T) 01252 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nao_orb, nao_orb,& 01253 1.0_dp,admm_env%A,admm_env%work_orb_orb,0.0_dp,& 01254 admm_env%work_aux_orb,error) 01255 ! *** A*C*Y^(T)*C^(T)*A^(T) add to aux_aux, minus sign cancels 01256 CALL cp_fm_gemm('N', 'T', nao_aux_fit, nao_aux_fit, nao_orb,& 01257 1.0_dp,admm_env%work_aux_orb,admm_env%A,1.0_dp,& 01258 admm_env%work_aux_aux,error) 01259 END IF 01260 01261 ! *** copy to sparse matrix 01262 CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, matrix_w_s,keep_sparsity=.TRUE.,& 01263 error=error) 01264 01265 01266 ! *** This can be done in one call w_total = w_alpha + w_beta 01267 CALL build_overlap_matrix(qs_env, nderivative=-1, & 01268 basis_set_id_a=use_aux_fit_basis_set, basis_set_id_b=use_aux_fit_basis_set, & 01269 sab_nl=qs_env%sab_aux_fit_asymm, calculate_forces=.TRUE., admm=.TRUE., & 01270 matrix_p=matrix_w_s, error=error) 01271 CALL build_overlap_matrix(qs_env, nderivative=-1, & 01272 basis_set_id_a=use_aux_fit_basis_set, basis_set_id_b=use_orb_basis_set, & 01273 sab_nl=qs_env%sab_aux_fit_vs_orb, calculate_forces=.TRUE., admm=.TRUE., & 01274 matrix_p=matrix_w_q, error=error) 01275 01276 01277 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 01278 qs_env%input,"DFT%PRINT%AO_MATRICES/W_MATRIX_AUX_FIT",error=error),cp_p_file)) THEN 01279 iw = cp_print_key_unit_nr(logger,qs_env%input,"DFT%PRINT%AO_MATRICES/W_MATRIX_AUX_FIT",& 01280 extension=".Log",error=error) 01281 CALL cp_dbcsr_write_sparse_matrix(matrix_w_s,4,6,qs_env,para_env,output_unit=iw,error=error) 01282 CALL cp_print_key_finished_output(iw,logger,qs_env%input,& 01283 "DFT%PRINT%AO_MATRICES/W_MATRIX_AUX_FIT", error=error) 01284 END IF 01285 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 01286 qs_env%input,"DFT%PRINT%AO_MATRICES/W_MATRIX_AUX_FIT",error=error),cp_p_file)) THEN 01287 iw = cp_print_key_unit_nr(logger,qs_env%input,"DFT%PRINT%AO_MATRICES/W_MATRIX_AUX_FIT",& 01288 extension=".Log",error=error) 01289 CALL cp_dbcsr_write_sparse_matrix(matrix_w_q,4,6,qs_env,para_env,output_unit=iw,error=error) 01290 CALL cp_print_key_finished_output(iw,logger,qs_env%input,& 01291 "DFT%PRINT%AO_MATRICES/W_MATRIX_AUX_FIT", error=error) 01292 END IF 01293 01294 ! *** Deallocated weighted density matrices 01295 CALL cp_dbcsr_deallocate_matrix(matrix_w_s,error) 01296 CALL cp_dbcsr_deallocate_matrix(matrix_w_q,error) 01297 01298 CALL timestop(handle) 01299 01300 END SUBROUTINE calc_mixed_overlap_force 01301 01302 ! ***************************************************************************** 01328 SUBROUTINE create_admm_xc_section(qs_env, xc_section, admm_env,error) 01329 TYPE(qs_environment_type), POINTER :: qs_env 01330 TYPE(section_vals_type), POINTER :: xc_section 01331 TYPE(admm_type), POINTER :: admm_env 01332 TYPE(cp_error_type), INTENT(INOUT) :: error 01333 01334 CHARACTER(LEN=*), PARAMETER :: routineN = 'create_admm_xc_section', 01335 routineP = moduleN//':'//routineN 01336 01337 INTEGER :: hfx_potential_type, ifun, nfun 01338 LOGICAL :: failure, funct_found 01339 REAL(dp) :: cutoff_radius, hfx_fraction, 01340 omega, scale_x 01341 TYPE(section_vals_type), POINTER :: xc_fun, xc_fun_section 01342 01343 NULLIFY(admm_env%xc_section_aux, admm_env%xc_section_primary) 01344 failure=.FALSE. 01345 01346 !! ** Duplicate existing xc-section 01347 01348 CALL section_vals_duplicate(xc_section,admm_env%xc_section_aux,error=error) 01349 CALL section_vals_duplicate(xc_section,admm_env%xc_section_primary,error=error) 01350 !** Now modify the auxiliary basis 01351 !** First remove all functionals 01352 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_aux,"XC_FUNCTIONAL",error=error) 01353 01354 !* Overwrite possible shortcut 01355 CALL section_vals_val_set(xc_fun_section,"_SECTION_PARAMETERS_",& 01356 i_val=xc_funct_no_shortcut,error=error) 01357 01358 !** Get number of Functionals in the list 01359 ifun = 0 01360 nfun = 0 01361 DO 01362 ifun = ifun+1 01363 xc_fun => section_vals_get_subs_vals2(xc_fun_section,i_section=ifun,error=error) 01364 IF (.NOT.ASSOCIATED(xc_fun)) EXIT 01365 nfun = nfun + 1 01366 END DO 01367 01368 ifun = 0 01369 DO ifun = 1,nfun 01370 xc_fun => section_vals_get_subs_vals2(xc_fun_section,i_section=1,error=error) 01371 IF (.NOT.ASSOCIATED(xc_fun)) EXIT 01372 CALL section_vals_remove_values(xc_fun, error=error) 01373 END DO 01374 01375 hfx_potential_type = qs_env%x_data(1,1)%potential_parameter%potential_type 01376 hfx_fraction = qs_env%x_data(1,1)%general_parameter%fraction 01377 01378 !! ** Add functionals evaluated with auxiliary basis 01379 SELECT CASE (hfx_potential_type) 01380 CASE (do_hfx_potential_coulomb) 01381 CALL section_vals_val_set(xc_fun_section,"PBE%_SECTION_PARAMETERS_",& 01382 l_val=.TRUE.,error=error) 01383 CALL section_vals_val_set(xc_fun_section,"PBE%SCALE_X",& 01384 r_val=-hfx_fraction,error=error) 01385 CALL section_vals_val_set(xc_fun_section,"PBE%SCALE_C",& 01386 r_val=0.0_dp,error=error) 01387 CASE (do_hfx_potential_short) 01388 omega = qs_env%x_data(1,1)%potential_parameter%omega 01389 CALL section_vals_val_set(xc_fun_section,"XWPBE%_SECTION_PARAMETERS_",& 01390 l_val=.TRUE.,error=error) 01391 CALL section_vals_val_set(xc_fun_section,"XWPBE%SCALE_X",& 01392 r_val=-hfx_fraction,error=error) 01393 CALL section_vals_val_set(xc_fun_section,"XWPBE%SCALE_X0",& 01394 r_val=0.0_dp,error=error) 01395 CALL section_vals_val_set(xc_fun_section,"XWPBE%OMEGA",& 01396 r_val=omega,error=error) 01397 CASE (do_hfx_potential_truncated) 01398 cutoff_radius = qs_env%x_data(1,1)%potential_parameter%cutoff_radius 01399 CALL section_vals_val_set(xc_fun_section,"PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_",& 01400 l_val=.TRUE.,error=error) 01401 CALL section_vals_val_set(xc_fun_section,"PBE_HOLE_T_C_LR%SCALE_X",& 01402 r_val=hfx_fraction,error=error) 01403 CALL section_vals_val_set(xc_fun_section,"PBE_HOLE_T_C_LR%CUTOFF_RADIUS",& 01404 r_val=cutoff_radius,error=error) 01405 CALL section_vals_val_set(xc_fun_section,"XWPBE%_SECTION_PARAMETERS_",& 01406 l_val=.TRUE.,error=error) 01407 CALL section_vals_val_set(xc_fun_section,"XWPBE%SCALE_X",& 01408 r_val=0.0_dp,error=error) 01409 CALL section_vals_val_set(xc_fun_section,"XWPBE%SCALE_X0",& 01410 r_val=-hfx_fraction,error=error) 01411 CASE DEFAULT 01412 CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) 01413 END SELECT 01414 01415 01416 !** Now modify the functionals for the primary basis 01417 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary,"XC_FUNCTIONAL",error=error) 01418 !* Overwrite possible shortcut 01419 CALL section_vals_val_set(xc_fun_section,"_SECTION_PARAMETERS_",& 01420 i_val=xc_funct_no_shortcut,error=error) 01421 01422 01423 SELECT CASE (hfx_potential_type) 01424 CASE (do_hfx_potential_coulomb) 01425 ifun = 0 01426 funct_found = .FALSE. 01427 DO 01428 ifun = ifun+1 01429 xc_fun => section_vals_get_subs_vals2(xc_fun_section,i_section=ifun,error=error) 01430 IF (.NOT.ASSOCIATED(xc_fun)) EXIT 01431 IF( xc_fun%section%name == "PBE" ) THEN 01432 funct_found = .TRUE. 01433 END IF 01434 END DO 01435 IF( .NOT. funct_found ) THEN 01436 CALL section_vals_val_set(xc_fun_section,"PBE%_SECTION_PARAMETERS_",& 01437 l_val=.TRUE.,error=error) 01438 CALL section_vals_val_set(xc_fun_section,"PBE%SCALE_X",& 01439 r_val=hfx_fraction,error=error) 01440 CALL section_vals_val_set(xc_fun_section,"PBE%SCALE_C",& 01441 r_val=0.0_dp,error=error) 01442 ELSE 01443 CALL section_vals_val_get(xc_fun_section,"PBE%SCALE_X",& 01444 r_val=scale_x,error=error) 01445 scale_x = scale_x + hfx_fraction 01446 CALL section_vals_val_set(xc_fun_section,"PBE%SCALE_X",& 01447 r_val=scale_x,error=error) 01448 END IF 01449 CASE (do_hfx_potential_short) 01450 omega = qs_env%x_data(1,1)%potential_parameter%omega 01451 ifun = 0 01452 funct_found = .FALSE. 01453 DO 01454 ifun = ifun+1 01455 xc_fun => section_vals_get_subs_vals2(xc_fun_section,i_section=ifun,error=error) 01456 IF (.NOT.ASSOCIATED(xc_fun)) EXIT 01457 IF( xc_fun%section%name == "XWPBE" ) THEN 01458 funct_found = .TRUE. 01459 END IF 01460 END DO 01461 IF( .NOT. funct_found ) THEN 01462 CALL section_vals_val_set(xc_fun_section,"XWPBE%_SECTION_PARAMETERS_",& 01463 l_val=.TRUE.,error=error) 01464 CALL section_vals_val_set(xc_fun_section,"XWPBE%SCALE_X",& 01465 r_val=hfx_fraction,error=error) 01466 CALL section_vals_val_set(xc_fun_section,"XWPBE%SCALE_X0",& 01467 r_val=0.0_dp,error=error) 01468 CALL section_vals_val_set(xc_fun_section,"XWPBE%OMEGA",& 01469 r_val=omega,error=error) 01470 ELSE 01471 CALL section_vals_val_get(xc_fun_section,"XWPBE%SCALE_X",& 01472 r_val=scale_x,error=error) 01473 scale_x = scale_x + hfx_fraction 01474 CALL section_vals_val_set(xc_fun_section,"XWPBE%SCALE_X",& 01475 r_val=scale_x,error=error) 01476 END IF 01477 CASE (do_hfx_potential_truncated) 01478 cutoff_radius = qs_env%x_data(1,1)%potential_parameter%cutoff_radius 01479 ifun = 0 01480 funct_found = .FALSE. 01481 DO 01482 ifun = ifun+1 01483 xc_fun => section_vals_get_subs_vals2(xc_fun_section,i_section=ifun,error=error) 01484 IF (.NOT.ASSOCIATED(xc_fun)) EXIT 01485 IF( xc_fun%section%name == "PBE_HOLE_T_C_LR" ) THEN 01486 funct_found = .TRUE. 01487 END IF 01488 END DO 01489 IF( .NOT. funct_found ) THEN 01490 CALL section_vals_val_set(xc_fun_section,"PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_",& 01491 l_val=.TRUE.,error=error) 01492 CALL section_vals_val_set(xc_fun_section,"PBE_HOLE_T_C_LR%SCALE_X",& 01493 r_val=-hfx_fraction,error=error) 01494 CALL section_vals_val_set(xc_fun_section,"PBE_HOLE_T_C_LR%CUTOFF_RADIUS",& 01495 r_val=cutoff_radius,error=error) 01496 01497 ELSE 01498 CALL section_vals_val_get(xc_fun_section,"PBE_HOLE_T_C_LR%SCALE_X",& 01499 r_val=scale_x,error=error) 01500 scale_x = scale_x - hfx_fraction 01501 CALL section_vals_val_set(xc_fun_section,"PBE_HOLE_T_C_LR%SCALE_X",& 01502 r_val=scale_x,error=error) 01503 CALL section_vals_val_set(xc_fun_section,"PBE_HOLE_T_C_LR%CUTOFF_RADIUS",& 01504 r_val=cutoff_radius,error=error) 01505 END IF 01506 ifun = 0 01507 funct_found = .FALSE. 01508 DO 01509 ifun = ifun+1 01510 xc_fun => section_vals_get_subs_vals2(xc_fun_section,i_section=ifun,error=error) 01511 IF (.NOT.ASSOCIATED(xc_fun)) EXIT 01512 IF( xc_fun%section%name == "XWPBE" ) THEN 01513 funct_found = .TRUE. 01514 END IF 01515 END DO 01516 IF( .NOT. funct_found ) THEN 01517 CALL section_vals_val_set(xc_fun_section,"XWPBE%_SECTION_PARAMETERS_",& 01518 l_val=.TRUE.,error=error) 01519 CALL section_vals_val_set(xc_fun_section,"XWPBE%SCALE_X0",& 01520 r_val=hfx_fraction,error=error) 01521 CALL section_vals_val_set(xc_fun_section,"XWPBE%SCALE_X",& 01522 r_val=0.0_dp,error=error) 01523 01524 ELSE 01525 CALL section_vals_val_get(xc_fun_section,"XWPBE%SCALE_X0",& 01526 r_val=scale_x,error=error) 01527 scale_x = scale_x + hfx_fraction 01528 CALL section_vals_val_set(xc_fun_section,"XWPBE%SCALE_X0",& 01529 r_val=scale_x,error=error) 01530 END IF 01531 01532 END SELECT 01533 01534 01535 IF( 1==0 ) THEN 01536 WRITE(*,*) "primary" 01537 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary,"XC_FUNCTIONAL",error=error) 01538 ifun = 0 01539 funct_found = .FALSE. 01540 DO 01541 ifun = ifun+1 01542 xc_fun => section_vals_get_subs_vals2(xc_fun_section,i_section=ifun,error=error) 01543 IF (.NOT.ASSOCIATED(xc_fun)) EXIT 01544 01545 scale_x=-1000.0_dp 01546 IF(xc_fun%section%name /= "LYP" .AND. xc_fun%section%name /= "VWN") THEN 01547 CALL section_vals_val_get(xc_fun,"SCALE_X",& 01548 r_val=scale_x,error=error) 01549 END IF 01550 IF( xc_fun%section%name == "XWPBE" ) THEN 01551 CALL section_vals_val_get(xc_fun,"SCALE_X0",& 01552 r_val=hfx_fraction,error=error) 01553 01554 WRITE(*,*) xc_fun%section%name, scale_x, hfx_fraction 01555 ELSE 01556 WRITE(*,*) xc_fun%section%name, scale_x 01557 END IF 01558 END DO 01559 01560 WRITE(*,*) "auxiliary" 01561 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_aux,"XC_FUNCTIONAL",error=error) 01562 ifun = 0 01563 funct_found = .FALSE. 01564 DO 01565 ifun = ifun+1 01566 xc_fun => section_vals_get_subs_vals2(xc_fun_section,i_section=ifun,error=error) 01567 IF (.NOT.ASSOCIATED(xc_fun)) EXIT 01568 scale_x=-1000.0_dp 01569 IF(xc_fun%section%name /= "LYP" .AND. xc_fun%section%name /= "VWN") THEN 01570 CALL section_vals_val_get(xc_fun,"SCALE_X",& 01571 r_val=scale_x,error=error) 01572 END IF 01573 IF( xc_fun%section%name == "XWPBE" ) THEN 01574 CALL section_vals_val_get(xc_fun,"SCALE_X0",& 01575 r_val=hfx_fraction,error=error) 01576 01577 WRITE(*,*) xc_fun%section%name, scale_x, hfx_fraction 01578 ELSE 01579 WRITE(*,*) xc_fun%section%name, scale_x 01580 END IF 01581 END DO 01582 END IF 01583 01584 END SUBROUTINE create_admm_xc_section 01585 01586 SUBROUTINE merge_ks_matrix_cauchy_subspace(ispin, admm_env, mo_set, mo_coeff, mo_coeff_aux_fit, & 01587 matrix_ks, matrix_ks_aux_fit, matrix_s, matrix_p_aux_fit, error) 01588 INTEGER, INTENT(IN) :: ispin 01589 TYPE(admm_type), POINTER :: admm_env 01590 TYPE(mo_set_type), POINTER :: mo_set 01591 TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit 01592 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01593 POINTER :: matrix_ks, matrix_ks_aux_fit, 01594 matrix_s, matrix_p_aux_fit 01595 TYPE(cp_error_type), INTENT(INOUT) :: error 01596 01597 CHARACTER(LEN=*), PARAMETER :: 01598 routineN = 'merge_ks_matrix_cauchy_subspace', 01599 routineP = moduleN//':'//routineN 01600 01601 INTEGER :: handle, nao_aux_fit, nao_orb, 01602 nmo 01603 INTEGER, SAVE :: counter = 0 01604 LOGICAL :: failure 01605 TYPE(cp_dbcsr_type), POINTER :: matrix_k_tilde 01606 01607 failure = .FALSE. 01608 01609 CALL timeset(routineN,handle) 01610 01611 counter = counter + 1 01612 nao_aux_fit = admm_env%nao_aux_fit 01613 nao_orb = admm_env%nao_orb 01614 nmo = admm_env%nmo(ispin) 01615 01616 !! Calculate Lambda^{-2} 01617 CALL cp_fm_to_fm(admm_env%lambda(ispin)%matrix, admm_env%work_nmo_nmo1(ispin)%matrix, error=error) 01618 CALL cp_fm_cholesky_decompose(admm_env%work_nmo_nmo1(ispin)%matrix,error=error) 01619 CALL cp_fm_cholesky_invert(admm_env%work_nmo_nmo1(ispin)%matrix,error=error) 01620 !! Symmetrize the guy 01621 CALL cp_fm_upper_to_full(admm_env%work_nmo_nmo1(ispin)%matrix,admm_env%lambda_inv2(ispin)%matrix,error=error) 01622 !! Take square 01623 CALL cp_fm_gemm('N', 'T', nmo, nmo, nmo,& 01624 1.0_dp,admm_env%work_nmo_nmo1(ispin)%matrix,admm_env%work_nmo_nmo1(ispin)%matrix,0.0_dp,& 01625 admm_env%lambda_inv2(ispin)%matrix,error) 01626 01627 !! ** C_hat = AC 01628 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nmo, nao_orb,& 01629 1.0_dp,admm_env%A,mo_coeff,0.0_dp,& 01630 admm_env%C_hat(ispin)%matrix,error) 01631 01632 !! calc P_tilde from C_hat 01633 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nmo, nmo,& 01634 1.0_dp,admm_env%C_hat(ispin)%matrix,admm_env%lambda_inv(ispin)%matrix,0.0_dp,& 01635 admm_env%work_aux_nmo(ispin)%matrix,error) 01636 01637 CALL cp_fm_gemm('N', 'T', nao_aux_fit, nao_aux_fit, nmo,& 01638 1.0_dp,admm_env%C_hat(ispin)%matrix, admm_env%work_aux_nmo(ispin)%matrix,0.0_dp,& 01639 admm_env%P_tilde(ispin)%matrix,error) 01640 01641 01642 !! ** C_hat*Lambda^{-2} 01643 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nmo, nmo,& 01644 1.0_dp,admm_env%C_hat(ispin)%matrix,admm_env%lambda_inv2(ispin)%matrix,0.0_dp,& 01645 admm_env%work_aux_nmo(ispin)%matrix,error) 01646 01647 !! ** C_hat*Lambda^{-2}*C_hat^T 01648 CALL cp_fm_gemm('N', 'T', nao_aux_fit, nao_aux_fit, nmo,& 01649 1.0_dp,admm_env%work_aux_nmo(ispin)%matrix,admm_env%C_hat(ispin)%matrix,0.0_dp,& 01650 admm_env%work_aux_aux,error) 01651 01652 01653 !! ** S*C_hat*Lambda^{-2}*C_hat^T 01654 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit,& 01655 1.0_dp,admm_env%S,admm_env%work_aux_aux,0.0_dp,& 01656 admm_env%work_aux_aux2,error) 01657 01658 01659 CALL copy_dbcsr_to_fm(matrix_ks_aux_fit(ispin)%matrix,admm_env%K(ispin)%matrix,error) 01660 CALL cp_fm_upper_to_full(admm_env%K(ispin)%matrix,admm_env%work_aux_aux,error=error) 01661 01662 !! ** S*C_hat*Lambda^{-2}*C_hat^T*H_tilde 01663 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit,& 01664 1.0_dp,admm_env%work_aux_aux2,admm_env%K(ispin)%matrix,0.0_dp,& 01665 admm_env%work_aux_aux,error) 01666 01667 !! ** P_tilde*S 01668 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit,& 01669 1.0_dp,admm_env%P_tilde(ispin)%matrix,admm_env%S,0.0_dp,& 01670 admm_env%work_aux_aux2,error) 01671 01672 01673 !! ** -S*C_hat*Lambda^{-2}*C_hat^T*H_tilde*P_tilde*S 01674 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit,& 01675 -1.0_dp,admm_env%work_aux_aux,admm_env%work_aux_aux2,0.0_dp,& 01676 admm_env%work_aux_aux3,error) 01677 01678 01679 !! ** -S*C_hat*Lambda^{-2}*C_hat^T*H_tilde*P_tilde*S+S*C_hat*Lambda^{-2}*C_hat^T*H_tilde 01680 CALL cp_fm_scale_and_add(1.0_dp,admm_env%work_aux_aux3,1.0_dp,admm_env%work_aux_aux,error) 01681 01682 01683 !! first_part*A 01684 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit,& 01685 1.0_dp,admm_env%work_aux_aux3,admm_env%A,0.0_dp,& 01686 admm_env%work_aux_orb,error) 01687 01688 !! + first_part^T*A 01689 CALL cp_fm_gemm('T', 'N', nao_aux_fit, nao_orb, nao_aux_fit,& 01690 1.0_dp,admm_env%work_aux_aux3,admm_env%A,1.0_dp,& 01691 admm_env%work_aux_orb,error) 01692 01693 01694 01695 !! A^T*(first+seccond)=H 01696 CALL cp_fm_gemm('T', 'N', nao_orb, nao_orb, nao_aux_fit,& 01697 1.0_dp,admm_env%A,admm_env%work_aux_orb,0.0_dp,& 01698 admm_env%work_orb_orb,error) 01699 01700 01701 NULLIFY(matrix_k_tilde) 01702 ALLOCATE(matrix_k_tilde) 01703 CALL cp_dbcsr_init (matrix_k_tilde, error) 01704 CALL cp_dbcsr_create(matrix_k_tilde, 'MATRIX K_tilde', & 01705 cp_dbcsr_distribution(matrix_ks(ispin)%matrix), dbcsr_type_symmetric, cp_dbcsr_row_block_sizes(matrix_ks(ispin)%matrix),& 01706 cp_dbcsr_col_block_sizes(matrix_ks(ispin)%matrix), & 01707 cp_dbcsr_get_num_blocks(matrix_ks(ispin)%matrix), & 01708 cp_dbcsr_get_data_size(matrix_ks(ispin)%matrix),& 01709 cp_dbcsr_get_data_type(matrix_ks(ispin)%matrix), error=error) 01710 01711 01712 CALL cp_fm_to_fm(admm_env%work_orb_orb, admm_env%ks_to_be_merged(ispin)%matrix, error=error) 01713 01714 CALL cp_dbcsr_copy(matrix_k_tilde, matrix_ks(ispin)%matrix, error=error) 01715 CALL cp_dbcsr_set(matrix_k_tilde, 0.0_dp, error) 01716 CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, matrix_k_tilde, keep_sparsity=.TRUE.,& 01717 error=error) 01718 01719 CALL cp_fm_gemm('N', 'N', nao_orb, nmo, nao_orb,& 01720 1.0_dp,admm_env%work_orb_orb,mo_coeff,0.0_dp,& 01721 admm_env%mo_derivs_tmp(ispin)%matrix,error) 01722 01723 01724 CALL cp_dbcsr_add(matrix_ks(ispin)%matrix, matrix_k_tilde, 1.0_dp, 1.0_dp, error) 01725 01726 CALL cp_dbcsr_deallocate_matrix(matrix_k_tilde,error) 01727 01728 CALL timestop(handle) 01729 01730 END SUBROUTINE merge_ks_matrix_cauchy_subspace 01731 01732 SUBROUTINE merge_ks_matrix_cauchy_blocked(ispin, admm_env, mo_set, mo_coeff, mo_coeff_aux_fit, & 01733 matrix_ks, matrix_ks_aux_fit, matrix_s, matrix_p_aux_fit,& 01734 matrix_p, error) 01735 INTEGER, INTENT(IN) :: ispin 01736 TYPE(admm_type), POINTER :: admm_env 01737 TYPE(mo_set_type), POINTER :: mo_set 01738 TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit 01739 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01740 POINTER :: matrix_ks, matrix_ks_aux_fit, 01741 matrix_s, matrix_p_aux_fit, 01742 matrix_p 01743 TYPE(cp_error_type), INTENT(INOUT) :: error 01744 01745 CHARACTER(LEN=*), PARAMETER :: 01746 routineN = 'merge_ks_matrix_cauchy_blocked', 01747 routineP = moduleN//':'//routineN 01748 01749 INTEGER :: blk, handle, i, iatom, j, 01750 jatom, nao_aux_fit, nao_orb, 01751 nmo, nspins 01752 INTEGER, SAVE :: counter = 0 01753 LOGICAL :: failure 01754 REAL(dp) :: eig_diff, pole, tmp 01755 REAL(dp), DIMENSION(:, :), POINTER :: sparse_block 01756 TYPE(cp_dbcsr_iterator) :: iter 01757 TYPE(cp_dbcsr_type), POINTER :: matrix_k_tilde 01758 01759 failure = .FALSE. 01760 01761 CALL timeset(routineN,handle) 01762 01763 counter = counter + 1 01764 nao_aux_fit = admm_env%nao_aux_fit 01765 nao_orb = admm_env%nao_orb 01766 nmo = admm_env%nmo(ispin) 01767 nspins = SIZE(admm_env%P_to_be_purified) 01768 01769 CALL cp_fm_to_fm(admm_env%S, admm_env%work_aux_aux, error=error) 01770 CALL cp_fm_to_fm(admm_env%P_to_be_purified(ispin)%matrix, admm_env%work_aux_aux2, error=error) 01771 01772 CALL cp_fm_cholesky_decompose(admm_env%work_aux_aux,error=error) 01773 01774 CALL cp_fm_cholesky_reduce(admm_env%work_aux_aux2, admm_env%work_aux_aux, itype=3, error=error) 01775 01776 CALL cp_fm_syevd(admm_env%work_aux_aux2,admm_env%R_purify(ispin)%matrix,& 01777 admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data,error=error) 01778 01779 CALL cp_fm_cholesky_restore(admm_env%R_purify(ispin)%matrix, nao_aux_fit,admm_env%work_aux_aux, & 01780 admm_env%work_aux_aux3,op="MULTIPLY",pos="LEFT", transa="T", error=error) 01781 01782 CALL cp_fm_to_fm(admm_env%work_aux_aux3, admm_env%R_purify(ispin)%matrix, error=error) 01783 01784 ! *** Construct Matrix M for Hadamard Product 01785 pole = 0.0_dp 01786 DO i=1,nao_aux_fit 01787 DO j=i,nao_aux_fit 01788 eig_diff = ( admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(i) -& 01789 admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(j) ) 01790 ! *** two eigenvalues could be the degenerated. In that case use 2nd order formula for the poles 01791 IF( ABS(eig_diff) == 0.0_dp ) THEN 01792 pole = delta(admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(i)-0.5_dp) 01793 CALL cp_fm_set_element(admm_env%M_purify(ispin)%matrix,i,j,pole,error) 01794 ELSE 01795 pole = 1.0_dp/(admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(i)-& 01796 admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(j)) 01797 tmp = Heaviside(admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(i)-0.5_dp) 01798 tmp = tmp - Heaviside(admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(j)-0.5_dp) 01799 pole = tmp*pole 01800 CALL cp_fm_set_element(admm_env%M_purify(ispin)%matrix,i,j,pole,error) 01801 END IF 01802 END DO 01803 END DO 01804 CALL cp_fm_upper_to_full(admm_env%M_purify(ispin)%matrix,admm_env%work_aux_aux,error=error) 01805 01806 CALL copy_dbcsr_to_fm(matrix_ks_aux_fit(ispin)%matrix,admm_env%K(ispin)%matrix,error) 01807 CALL cp_fm_upper_to_full(admm_env%K(ispin)%matrix,admm_env%work_aux_aux,error=error) 01808 01809 01810 !! S^(-1)*R 01811 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit,& 01812 1.0_dp,admm_env%S_inv,admm_env%R_purify(ispin)%matrix,0.0_dp,& 01813 admm_env%work_aux_aux,error) 01814 !! K*S^(-1)*R 01815 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit,& 01816 1.0_dp,admm_env%K(ispin)%matrix,admm_env%work_aux_aux,0.0_dp,& 01817 admm_env%work_aux_aux2,error) 01818 !! R^T*S^(-1)*K*S^(-1)*R 01819 CALL cp_fm_gemm('T', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit,& 01820 1.0_dp,admm_env%work_aux_aux,admm_env%work_aux_aux2,0.0_dp,& 01821 admm_env%work_aux_aux3,error) 01822 !! R^T*S^(-1)*K*S^(-1)*R x M 01823 CALL cp_fm_schur_product(admm_env%work_aux_aux3, admm_env%M_purify(ispin)%matrix,& 01824 admm_env%work_aux_aux,error) 01825 01826 !! R^T*A 01827 CALL cp_fm_gemm('T', 'N', nao_aux_fit, nao_orb, nao_aux_fit,& 01828 1.0_dp, admm_env%R_purify(ispin)%matrix, admm_env%A, 0.0_dp,& 01829 admm_env%work_aux_orb,error) 01830 01831 !! (R^T*S^(-1)*K*S^(-1)*R x M) * R^T*A 01832 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit,& 01833 1.0_dp, admm_env%work_aux_aux, admm_env%work_aux_orb, 0.0_dp,& 01834 admm_env%work_aux_orb2,error) 01835 !! A^T*R*(R^T*S^(-1)*K*S^(-1)*R x M) * R^T*A 01836 CALL cp_fm_gemm('T', 'N', nao_orb, nao_orb, nao_aux_fit,& 01837 1.0_dp, admm_env%work_aux_orb, admm_env%work_aux_orb2, 0.0_dp,& 01838 admm_env%work_orb_orb,error) 01839 01840 01841 NULLIFY(matrix_k_tilde) 01842 ALLOCATE(matrix_k_tilde) 01843 CALL cp_dbcsr_init (matrix_k_tilde, error) 01844 CALL cp_dbcsr_create(matrix_k_tilde, 'MATRIX K_tilde', & 01845 cp_dbcsr_distribution(matrix_ks(ispin)%matrix), dbcsr_type_symmetric,& 01846 cp_dbcsr_row_block_sizes(matrix_ks(ispin)%matrix),& 01847 cp_dbcsr_col_block_sizes(matrix_ks(ispin)%matrix), & 01848 cp_dbcsr_get_num_blocks(matrix_ks(ispin)%matrix), cp_dbcsr_get_data_size(matrix_ks(ispin)%matrix),& 01849 cp_dbcsr_get_data_type(matrix_ks(ispin)%matrix), & 01850 error=error) 01851 01852 CALL cp_fm_to_fm(admm_env%work_orb_orb, admm_env%ks_to_be_merged(ispin)%matrix, error=error) 01853 01854 CALL cp_dbcsr_copy(matrix_k_tilde, matrix_ks(ispin)%matrix, error=error) 01855 CALL cp_dbcsr_set(matrix_k_tilde, 0.0_dp, error) 01856 CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, matrix_k_tilde, keep_sparsity=.TRUE., & 01857 error=error) 01858 01859 ! ** now loop through the list and nullify blocks 01860 CALL cp_dbcsr_iterator_start(iter, matrix_k_tilde) 01861 DO WHILE (cp_dbcsr_iterator_blocks_left (iter)) 01862 CALL cp_dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block, blk) 01863 IF( admm_env%block_map(iatom,jatom) == 0 ) THEN 01864 sparse_block = 0.0_dp 01865 END IF 01866 END DO 01867 CALL cp_dbcsr_iterator_stop(iter) 01868 01869 CALL cp_dbcsr_add(matrix_ks(ispin)%matrix, matrix_k_tilde, 1.0_dp, 1.0_dp, error) 01870 01871 CALL cp_dbcsr_deallocate_matrix(matrix_k_tilde,error) 01872 01873 CALL timestop(handle) 01874 01875 END SUBROUTINE merge_ks_matrix_cauchy_blocked 01876 01877 01878 SUBROUTINE merge_ks_matrix_none(ispin, admm_env, mo_set, mo_coeff, mo_coeff_aux_fit, & 01879 matrix_ks, matrix_ks_aux_fit, matrix_s, matrix_p_aux_fit, error) 01880 INTEGER, INTENT(IN) :: ispin 01881 TYPE(admm_type), POINTER :: admm_env 01882 TYPE(mo_set_type), POINTER :: mo_set 01883 TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit 01884 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01885 POINTER :: matrix_ks, matrix_ks_aux_fit, 01886 matrix_s, matrix_p_aux_fit 01887 TYPE(cp_error_type), INTENT(INOUT) :: error 01888 01889 CHARACTER(LEN=*), PARAMETER :: routineN = 'merge_ks_matrix_none', 01890 routineP = moduleN//':'//routineN 01891 01892 INTEGER :: handle, nao_aux_fit, nao_orb, 01893 nmo 01894 INTEGER, SAVE :: counter = 0 01895 LOGICAL :: failure 01896 TYPE(cp_dbcsr_type), POINTER :: matrix_k_tilde 01897 01898 failure = .FALSE. 01899 01900 CALL timeset(routineN,handle) 01901 01902 counter = counter + 1 01903 nao_aux_fit = admm_env%nao_aux_fit 01904 nao_orb = admm_env%nao_orb 01905 nmo = admm_env%nmo(ispin) 01906 01907 CALL copy_dbcsr_to_fm(matrix_ks_aux_fit(ispin)%matrix,admm_env%K(ispin)%matrix,error) 01908 CALL cp_fm_upper_to_full(admm_env%K(ispin)%matrix,admm_env%work_aux_aux,error=error) 01909 01910 !! K*A 01911 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit,& 01912 1.0_dp,admm_env%K(ispin)%matrix,admm_env%A,0.0_dp,& 01913 admm_env%work_aux_orb,error) 01914 !! A^T*K*A 01915 CALL cp_fm_gemm('T', 'N', nao_orb, nao_orb, nao_aux_fit,& 01916 1.0_dp,admm_env%A,admm_env%work_aux_orb,0.0_dp,& 01917 admm_env%work_orb_orb,error) 01918 01919 01920 NULLIFY(matrix_k_tilde) 01921 ALLOCATE(matrix_k_tilde) 01922 CALL cp_dbcsr_init (matrix_k_tilde, error) 01923 CALL cp_dbcsr_create(matrix_k_tilde, 'MATRIX K_tilde', & 01924 cp_dbcsr_distribution(matrix_ks(ispin)%matrix), dbcsr_type_symmetric, cp_dbcsr_row_block_sizes(matrix_ks(ispin)%matrix),& 01925 cp_dbcsr_col_block_sizes(matrix_ks(ispin)%matrix), & 01926 cp_dbcsr_get_num_blocks(matrix_ks(ispin)%matrix), cp_dbcsr_get_data_size(matrix_ks(ispin)%matrix),& 01927 cp_dbcsr_get_data_type(matrix_ks(ispin)%matrix), error=error) 01928 01929 CALL cp_dbcsr_copy(matrix_k_tilde, matrix_ks(ispin)%matrix, error=error) 01930 CALL cp_dbcsr_set(matrix_k_tilde, 0.0_dp, error) 01931 CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, matrix_k_tilde, keep_sparsity=.TRUE.,& 01932 error=error) 01933 01934 CALL cp_dbcsr_add(matrix_ks(ispin)%matrix, matrix_k_tilde, 1.0_dp, 1.0_dp, error) 01935 01936 CALL cp_dbcsr_deallocate_matrix(matrix_k_tilde,error) 01937 01938 CALL timestop(handle) 01939 01940 END SUBROUTINE merge_ks_matrix_none 01941 01942 01943 SUBROUTINE print_matlab_matrix(matrix, filename, counter) 01944 TYPE(cp_fm_type), POINTER :: matrix 01945 CHARACTER(LEN=*) :: filename 01946 INTEGER :: counter 01947 01948 CHARACTER(LEN=20) :: string 01949 CHARACTER(LEN=50) :: filename_count 01950 INTEGER :: i, j 01951 01952 WRITE(filename_count,FMT='(A,I0)') filename, counter 01953 OPEN (unit=120,file=TRIM(filename_count)) 01954 j=SIZE(matrix%local_data,2) 01955 WRITE(string,FMT='(A1,I4,A7)') "(",j,"F20.16)" 01956 DO i=1,SIZE(matrix%local_data,1) 01957 WRITE(unit=120,FMT=string) matrix%local_data(i,:) 01958 END DO 01959 CLOSE (unit=120) 01960 END SUBROUTINE 01961 01962 SUBROUTINE remove_ks_matrix(ispin, admm_env, mo_set, mo_coeff, mo_coeff_aux_fit, & 01963 matrix_ks, matrix_ks_aux_fit, error) 01964 INTEGER, INTENT(IN) :: ispin 01965 TYPE(admm_type), POINTER :: admm_env 01966 TYPE(mo_set_type), POINTER :: mo_set 01967 TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit 01968 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01969 POINTER :: matrix_ks, matrix_ks_aux_fit 01970 TYPE(cp_error_type), INTENT(INOUT) :: error 01971 01972 CHARACTER(LEN=*), PARAMETER :: routineN = 'remove_ks_matrix', 01973 routineP = moduleN//':'//routineN 01974 01975 INTEGER :: handle 01976 LOGICAL :: failure 01977 TYPE(cp_dbcsr_type), POINTER :: matrix_k_tilde 01978 01979 failure = .FALSE. 01980 01981 CALL timeset(routineN,handle) 01982 01983 01984 NULLIFY(matrix_k_tilde) 01985 ALLOCATE(matrix_k_tilde) 01986 CALL cp_dbcsr_init (matrix_k_tilde, error) 01987 CALL cp_dbcsr_create(matrix_k_tilde, 'MATRIX K_tilde', & 01988 cp_dbcsr_distribution(matrix_ks(ispin)%matrix), dbcsr_type_symmetric, cp_dbcsr_row_block_sizes(matrix_ks(ispin)%matrix),& 01989 cp_dbcsr_col_block_sizes(matrix_ks(ispin)%matrix), cp_dbcsr_get_num_blocks(matrix_ks(ispin)%matrix), & 01990 cp_dbcsr_get_data_size(matrix_ks(ispin)%matrix),& 01991 cp_dbcsr_get_data_type(matrix_ks(ispin)%matrix), error=error) 01992 01993 01994 CALL cp_fm_to_fm(admm_env%work_orb_orb, admm_env%ks_to_be_merged(ispin)%matrix, error=error) 01995 01996 CALL cp_dbcsr_copy(matrix_k_tilde, matrix_ks(ispin)%matrix, error=error) 01997 CALL cp_dbcsr_set(matrix_k_tilde, 0.0_dp, error) 01998 CALL copy_fm_to_dbcsr(admm_env%ks_to_be_merged(ispin)%matrix, matrix_k_tilde, keep_sparsity=.TRUE.,& 01999 error=error) 02000 02001 CALL cp_dbcsr_add(matrix_ks(ispin)%matrix, matrix_k_tilde, 1.0_dp, -1.0_dp, error) 02002 02003 CALL cp_dbcsr_deallocate_matrix(matrix_k_tilde,error) 02004 02005 CALL timestop(handle) 02006 02007 END SUBROUTINE remove_ks_matrix 02008 02009 SUBROUTINE calculate_dm_mo_no_diag(admm_env,mo_set,density_matrix,ispin,error) 02010 TYPE(admm_type), POINTER :: admm_env 02011 TYPE(mo_set_type), POINTER :: mo_set 02012 TYPE(cp_dbcsr_type), POINTER :: density_matrix 02013 INTEGER :: ispin 02014 TYPE(cp_error_type), INTENT(inout) :: error 02015 02016 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_dm_mo_no_diag', 02017 routineP = moduleN//':'//routineN 02018 02019 INTEGER :: handle, nao_aux_fit, nmo 02020 REAL(KIND=dp) :: alpha 02021 02022 CALL timeset(routineN,handle) 02023 02024 02025 CALL cp_dbcsr_set(density_matrix,0.0_dp,error=error) 02026 nao_aux_fit = admm_env%nao_aux_fit 02027 nmo = admm_env%nmo(ispin) 02028 CALL cp_fm_to_fm(admm_env%C_hat(ispin)%matrix, admm_env%work_aux_nmo(ispin)%matrix, error=error) 02029 CALL cp_fm_column_scale(admm_env%work_aux_nmo(ispin)%matrix,mo_set%occupation_numbers(1:mo_set%homo)) 02030 02031 CALL cp_fm_gemm('N','N',nao_aux_fit,nmo,nmo,& 02032 1.0_dp,admm_env%work_aux_nmo(ispin)%matrix,admm_env%lambda_inv(ispin)%matrix,0.0_dp,& 02033 admm_env%work_aux_nmo2(ispin)%matrix,error) 02034 02035 02036 IF ( .NOT. mo_set%uniform_occupation ) THEN ! not all orbitals 1..homo are equally occupied 02037 alpha=1.0_dp 02038 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix,& 02039 matrix_v=admm_env%C_hat(ispin)%matrix,& 02040 matrix_g=admm_env%work_aux_nmo2(ispin)%matrix,& 02041 ncol=mo_set%homo,& 02042 alpha=alpha,error=error) 02043 ELSE 02044 alpha=1.0_dp 02045 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix,& 02046 matrix_v=admm_env%C_hat(ispin)%matrix,& 02047 matrix_g=admm_env%work_aux_nmo2(ispin)%matrix,& 02048 ncol=mo_set%homo,& 02049 alpha=alpha,error=error) 02050 ENDIF 02051 02052 CALL timestop(handle) 02053 02054 END SUBROUTINE calculate_dm_mo_no_diag 02055 02056 SUBROUTINE purify_dm_cauchy_blocked(admm_env,mo_set,density_matrix,ispin,error) 02057 02058 TYPE(admm_type), POINTER :: admm_env 02059 TYPE(mo_set_type), POINTER :: mo_set 02060 TYPE(cp_dbcsr_type), POINTER :: density_matrix 02061 INTEGER :: ispin 02062 TYPE(cp_error_type), INTENT(inout) :: error 02063 02064 CHARACTER(len=*), PARAMETER :: routineN = 'purify_dm_cauchy_blocked', 02065 routineP = moduleN//':'//routineN 02066 02067 INTEGER :: handle, i, nao_aux_fit, 02068 nao_orb, nmo, nspins 02069 REAL(KIND=dp) :: pole 02070 TYPE(cp_fm_type), POINTER :: mo_coeff_aux_fit 02071 02072 CALL timeset(routineN,handle) 02073 02074 nao_aux_fit = admm_env%nao_aux_fit 02075 nao_orb = admm_env%nao_orb 02076 nmo = admm_env%nmo(ispin) 02077 02078 nspins = SIZE(admm_env%P_to_be_purified) 02079 02080 CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff_aux_fit) 02081 02082 !! * For the time beeing, get the P to be purified from the mo_coeffs 02083 !! * This needs to be replaced with the a block modified P 02084 02085 ! CALL cp_fm_gemm('N','T',nao_aux_fit,nao_aux_fit,nmo,& 02086 ! 1.0_dp,mo_coeff_aux_fit,mo_coeff_aux_fit,0.0_dp,& 02087 ! admm_env%P_to_be_purified(ispin)%matrix,error) 02088 02089 CALL cp_fm_to_fm(admm_env%S, admm_env%work_aux_aux, error=error) 02090 CALL cp_fm_to_fm(admm_env%P_to_be_purified(ispin)%matrix, admm_env%work_aux_aux2, error=error) 02091 02092 CALL cp_fm_cholesky_decompose(admm_env%work_aux_aux,error=error) 02093 02094 CALL cp_fm_cholesky_reduce(admm_env%work_aux_aux2, admm_env%work_aux_aux, itype=3, error=error) 02095 02096 CALL cp_fm_syevd(admm_env%work_aux_aux2,admm_env%R_purify(ispin)%matrix,& 02097 admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data,error=error) 02098 02099 CALL cp_fm_cholesky_restore(admm_env%R_purify(ispin)%matrix, nao_aux_fit,admm_env%work_aux_aux, & 02100 admm_env%work_aux_aux3,op="MULTIPLY",pos="LEFT", transa="T", error=error) 02101 02102 CALL cp_fm_to_fm(admm_env%work_aux_aux3, admm_env%R_purify(ispin)%matrix, error=error) 02103 02104 ! *** Construct Matrix M for Hadamard Product 02105 CALL cp_fm_set_all(admm_env%M_purify(ispin)%matrix,0.0_dp,error=error) 02106 pole = 0.0_dp 02107 DO i=1,nao_aux_fit 02108 pole = Heaviside(admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(i)-0.5_dp) 02109 CALL cp_fm_set_element(admm_env%M_purify(ispin)%matrix,i,i,pole,error) 02110 END DO 02111 CALL cp_fm_upper_to_full(admm_env%M_purify(ispin)%matrix,admm_env%work_aux_aux,error=error) 02112 02113 CALL copy_dbcsr_to_fm(density_matrix,admm_env%work_aux_aux3,error) 02114 CALL cp_fm_upper_to_full(admm_env%work_aux_aux3,admm_env%work_aux_aux,error=error) 02115 02116 ! ** S^(-1)*R 02117 CALL cp_fm_gemm('N','N',nao_aux_fit,nao_aux_fit,nao_aux_fit,& 02118 1.0_dp,admm_env%S_inv,admm_env%R_purify(ispin)%matrix,0.0_dp,& 02119 admm_env%work_aux_aux,error) 02120 ! ** S^(-1)*R*M 02121 CALL cp_fm_gemm('N','N',nao_aux_fit,nao_aux_fit,nao_aux_fit,& 02122 1.0_dp,admm_env%work_aux_aux,admm_env%M_purify(ispin)%matrix,0.0_dp,& 02123 admm_env%work_aux_aux2,error) 02124 ! ** S^(-1)*R*M*R^T*S^(-1) 02125 CALL cp_fm_gemm('N','T',nao_aux_fit,nao_aux_fit,nao_aux_fit,& 02126 1.0_dp,admm_env%work_aux_aux2,admm_env%work_aux_aux,0.0_dp,& 02127 admm_env%work_aux_aux3,error) 02128 02129 CALL copy_fm_to_dbcsr(admm_env%work_aux_aux3, density_matrix,keep_sparsity=.TRUE., error=error) 02130 02131 IF( nspins == 1 ) THEN 02132 CALL cp_dbcsr_scale(density_matrix, 2.0_dp, error=error) 02133 END IF 02134 02135 CALL timestop(handle) 02136 02137 END SUBROUTINE purify_dm_cauchy_blocked 02138 02139 SUBROUTINE blockify_density_matrix(admm_env,density_matrix, density_matrix_aux,& 02140 ispin, nspins, error) 02141 TYPE(admm_type), POINTER :: admm_env 02142 TYPE(cp_dbcsr_type), POINTER :: density_matrix, 02143 density_matrix_aux 02144 INTEGER :: ispin, nspins 02145 TYPE(cp_error_type), INTENT(inout) :: error 02146 02147 CHARACTER(len=*), PARAMETER :: routineN = 'blockify_density_matrix', 02148 routineP = moduleN//':'//routineN 02149 02150 INTEGER :: blk, handle, iatom, jatom 02151 LOGICAL :: failure, found 02152 REAL(dp), DIMENSION(:, :), POINTER :: sparse_block, sparse_block_aux 02153 TYPE(cp_dbcsr_iterator) :: iter 02154 02155 failure = .FALSE. 02156 02157 CALL timeset(routineN,handle) 02158 02159 ! ** set blocked density matrix to 0 02160 CALL cp_dbcsr_set(density_matrix_aux, 0.0_dp, error) 02161 02162 ! ** now loop through the list and copy corresponding blocks 02163 CALL cp_dbcsr_iterator_start(iter, density_matrix) 02164 DO WHILE (cp_dbcsr_iterator_blocks_left (iter)) 02165 CALL cp_dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block, blk) 02166 IF( admm_env%block_map(iatom,jatom) == 1 ) THEN 02167 CALL cp_dbcsr_get_block_p(density_matrix_aux,& 02168 row=iatom,col=jatom,BLOCK=sparse_block_aux,found=found) 02169 IF( found ) THEN 02170 sparse_block_aux = sparse_block 02171 END IF 02172 02173 END IF 02174 END DO 02175 CALL cp_dbcsr_iterator_stop(iter) 02176 02177 CALL copy_dbcsr_to_fm(density_matrix_aux,admm_env%P_to_be_purified(ispin)%matrix,error) 02178 CALL cp_fm_upper_to_full(admm_env%P_to_be_purified(ispin)%matrix, admm_env%work_orb_orb2,error=error) 02179 02180 IF( nspins == 1 ) THEN 02181 CALL cp_fm_scale(0.5_dp, admm_env%P_to_be_purified(ispin)%matrix, error) 02182 END IF 02183 02184 CALL timestop(handle) 02185 END SUBROUTINE blockify_density_matrix 02186 02187 SUBROUTINE merge_ks_matrix_cauchy(ispin, admm_env, mo_set, mo_coeff, mo_coeff_aux_fit, & 02188 matrix_ks, matrix_ks_aux_fit, matrix_s, matrix_p_aux_fit, matrix_p, error) 02189 INTEGER, INTENT(IN) :: ispin 02190 TYPE(admm_type), POINTER :: admm_env 02191 TYPE(mo_set_type), POINTER :: mo_set 02192 TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit 02193 TYPE(cp_dbcsr_p_type), DIMENSION(:), 02194 POINTER :: matrix_ks, matrix_ks_aux_fit, 02195 matrix_s, matrix_p_aux_fit, 02196 matrix_p 02197 TYPE(cp_error_type), INTENT(INOUT) :: error 02198 02199 CHARACTER(LEN=*), PARAMETER :: routineN = 'merge_ks_matrix_cauchy', 02200 routineP = moduleN//':'//routineN 02201 02202 INTEGER :: handle, i, j, nao_aux_fit, 02203 nao_orb, nmo 02204 INTEGER, SAVE :: counter = 0 02205 LOGICAL :: failure 02206 REAL(dp) :: eig_diff, pole, tmp 02207 TYPE(cp_dbcsr_type), POINTER :: matrix_k_tilde 02208 02209 failure = .FALSE. 02210 02211 CALL timeset(routineN,handle) 02212 02213 counter = counter + 1 02214 nao_aux_fit = admm_env%nao_aux_fit 02215 nao_orb = admm_env%nao_orb 02216 nmo = admm_env%nmo(ispin) 02217 02218 !** Get P from mo_coeffs, otherwise we have troubles with occupation numbers ... 02219 CALL cp_fm_gemm('N', 'T', nao_orb, nao_orb, nmo,& 02220 1.0_dp, mo_coeff, mo_coeff, 0.0_dp,& 02221 admm_env%work_orb_orb,error) 02222 02223 !! A*P 02224 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nao_orb, nao_orb,& 02225 1.0_dp, admm_env%A, admm_env%work_orb_orb, 0.0_dp,& 02226 admm_env%work_aux_orb2,error) 02227 !! A*P*A^T 02228 CALL cp_fm_gemm('N', 'T', nao_aux_fit, nao_aux_fit, nao_orb,& 02229 1.0_dp, admm_env%work_aux_orb2, admm_env%A, 0.0_dp,& 02230 admm_env%P_to_be_purified(ispin)%matrix,error) 02231 02232 02233 CALL cp_fm_to_fm(admm_env%S, admm_env%work_aux_aux, error=error) 02234 CALL cp_fm_to_fm(admm_env%P_to_be_purified(ispin)%matrix, admm_env%work_aux_aux2, error=error) 02235 02236 CALL cp_fm_cholesky_decompose(admm_env%work_aux_aux,error=error) 02237 02238 CALL cp_fm_cholesky_reduce(admm_env%work_aux_aux2, admm_env%work_aux_aux, itype=3, error=error) 02239 02240 CALL cp_fm_syevd(admm_env%work_aux_aux2,admm_env%R_purify(ispin)%matrix,& 02241 admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data,error=error) 02242 02243 CALL cp_fm_cholesky_restore(admm_env%R_purify(ispin)%matrix, nao_aux_fit,admm_env%work_aux_aux, & 02244 admm_env%work_aux_aux3,op="MULTIPLY",pos="LEFT", transa="T", error=error) 02245 02246 CALL cp_fm_to_fm(admm_env%work_aux_aux3, admm_env%R_purify(ispin)%matrix, error=error) 02247 02248 ! *** Construct Matrix M for Hadamard Product 02249 pole = 0.0_dp 02250 DO i=1,nao_aux_fit 02251 DO j=i,nao_aux_fit 02252 eig_diff = ( admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(i) -& 02253 admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(j) ) 02254 ! *** two eigenvalues could be the degenerated. In that case use 2nd order formula for the poles 02255 IF( ABS(eig_diff) == 0.0_dp ) THEN 02256 pole = delta(admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(i)-0.5_dp) 02257 CALL cp_fm_set_element(admm_env%M_purify(ispin)%matrix,i,j,pole,error) 02258 ELSE 02259 pole = 1.0_dp/(admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(i)-& 02260 admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(j)) 02261 tmp = Heaviside(admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(i)-0.5_dp) 02262 tmp = tmp - Heaviside(admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(j)-0.5_dp) 02263 pole = tmp*pole 02264 CALL cp_fm_set_element(admm_env%M_purify(ispin)%matrix,i,j,pole,error) 02265 END IF 02266 END DO 02267 END DO 02268 CALL cp_fm_upper_to_full(admm_env%M_purify(ispin)%matrix,admm_env%work_aux_aux,error=error) 02269 02270 CALL copy_dbcsr_to_fm(matrix_ks_aux_fit(ispin)%matrix,admm_env%K(ispin)%matrix,error) 02271 CALL cp_fm_upper_to_full(admm_env%K(ispin)%matrix,admm_env%work_aux_aux,error=error) 02272 02273 02274 !! S^(-1)*R 02275 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit,& 02276 1.0_dp,admm_env%S_inv,admm_env%R_purify(ispin)%matrix,0.0_dp,& 02277 admm_env%work_aux_aux,error) 02278 !! K*S^(-1)*R 02279 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit,& 02280 1.0_dp,admm_env%K(ispin)%matrix,admm_env%work_aux_aux,0.0_dp,& 02281 admm_env%work_aux_aux2,error) 02282 !! R^T*S^(-1)*K*S^(-1)*R 02283 CALL cp_fm_gemm('T', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit,& 02284 1.0_dp,admm_env%work_aux_aux,admm_env%work_aux_aux2,0.0_dp,& 02285 admm_env%work_aux_aux3,error) 02286 !! R^T*S^(-1)*K*S^(-1)*R x M 02287 CALL cp_fm_schur_product(admm_env%work_aux_aux3, admm_env%M_purify(ispin)%matrix,& 02288 admm_env%work_aux_aux,error) 02289 02290 !! R^T*A 02291 CALL cp_fm_gemm('T', 'N', nao_aux_fit, nao_orb, nao_aux_fit,& 02292 1.0_dp, admm_env%R_purify(ispin)%matrix, admm_env%A, 0.0_dp,& 02293 admm_env%work_aux_orb,error) 02294 02295 !! (R^T*S^(-1)*K*S^(-1)*R x M) * R^T*A 02296 CALL cp_fm_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit,& 02297 1.0_dp, admm_env%work_aux_aux, admm_env%work_aux_orb, 0.0_dp,& 02298 admm_env%work_aux_orb2,error) 02299 !! A^T*R*(R^T*S^(-1)*K*S^(-1)*R x M) * R^T*A 02300 CALL cp_fm_gemm('T', 'N', nao_orb, nao_orb, nao_aux_fit,& 02301 1.0_dp, admm_env%work_aux_orb, admm_env%work_aux_orb2, 0.0_dp,& 02302 admm_env%work_orb_orb,error) 02303 02304 02305 NULLIFY(matrix_k_tilde) 02306 ALLOCATE(matrix_k_tilde) 02307 CALL cp_dbcsr_init (matrix_k_tilde, error) 02308 CALL cp_dbcsr_create(matrix_k_tilde, 'MATRIX K_tilde', & 02309 cp_dbcsr_distribution(matrix_ks(ispin)%matrix), dbcsr_type_symmetric, cp_dbcsr_row_block_sizes(matrix_ks(ispin)%matrix),& 02310 cp_dbcsr_col_block_sizes(matrix_ks(ispin)%matrix), cp_dbcsr_get_num_blocks(matrix_ks(ispin)%matrix),& 02311 cp_dbcsr_get_data_size( matrix_ks(ispin)%matrix),& 02312 cp_dbcsr_get_data_type(matrix_ks(ispin)%matrix), error=error) 02313 02314 CALL cp_fm_to_fm(admm_env%work_orb_orb, admm_env%ks_to_be_merged(ispin)%matrix, error=error) 02315 02316 CALL cp_dbcsr_copy(matrix_k_tilde, matrix_ks(ispin)%matrix, error=error) 02317 CALL cp_dbcsr_set(matrix_k_tilde, 0.0_dp, error) 02318 CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, matrix_k_tilde, keep_sparsity=.TRUE., error=error) 02319 02320 CALL cp_dbcsr_add(matrix_ks(ispin)%matrix, matrix_k_tilde, 1.0_dp, 1.0_dp, error) 02321 02322 CALL cp_dbcsr_deallocate_matrix(matrix_k_tilde,error) 02323 02324 CALL timestop(handle) 02325 02326 END SUBROUTINE merge_ks_matrix_cauchy 02327 02328 SUBROUTINE purify_density_matrix_cauchy(admm_env,mo_set,density_matrix,ispin,error) 02329 02330 TYPE(admm_type), POINTER :: admm_env 02331 TYPE(mo_set_type), POINTER :: mo_set 02332 TYPE(cp_dbcsr_type), POINTER :: density_matrix 02333 INTEGER :: ispin 02334 TYPE(cp_error_type), INTENT(inout) :: error 02335 02336 CHARACTER(len=*), PARAMETER :: routineN = 'purify_density_matrix_cauchy', 02337 routineP = moduleN//':'//routineN 02338 02339 INTEGER :: handle, i, nao_aux_fit, 02340 nao_orb, nmo, nspins 02341 REAL(KIND=dp) :: pole 02342 TYPE(cp_fm_type), POINTER :: mo_coeff_aux_fit 02343 02344 CALL timeset(routineN,handle) 02345 02346 nao_aux_fit = admm_env%nao_aux_fit 02347 nao_orb = admm_env%nao_orb 02348 nmo = admm_env%nmo(ispin) 02349 02350 nspins = SIZE(admm_env%P_to_be_purified) 02351 02352 CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff_aux_fit) 02353 02354 !! * For the time beeing, get the P to be purified from the mo_coeffs 02355 !! * This needs to be replaced with the a block modified P 02356 02357 CALL cp_fm_gemm('N','T',nao_aux_fit,nao_aux_fit,nmo,& 02358 1.0_dp,mo_coeff_aux_fit,mo_coeff_aux_fit,0.0_dp,& 02359 admm_env%P_to_be_purified(ispin)%matrix,error) 02360 02361 CALL cp_fm_to_fm(admm_env%S, admm_env%work_aux_aux, error=error) 02362 CALL cp_fm_to_fm(admm_env%P_to_be_purified(ispin)%matrix, admm_env%work_aux_aux2, error=error) 02363 02364 CALL cp_fm_cholesky_decompose(admm_env%work_aux_aux,error=error) 02365 02366 CALL cp_fm_cholesky_reduce(admm_env%work_aux_aux2, admm_env%work_aux_aux, itype=3, error=error) 02367 02368 CALL cp_fm_syevd(admm_env%work_aux_aux2,admm_env%R_purify(ispin)%matrix,& 02369 admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data,error=error) 02370 02371 CALL cp_fm_cholesky_restore(admm_env%R_purify(ispin)%matrix, nao_aux_fit,admm_env%work_aux_aux, & 02372 admm_env%work_aux_aux3,op="MULTIPLY",pos="LEFT", transa="T", error=error) 02373 02374 CALL cp_fm_to_fm(admm_env%work_aux_aux3, admm_env%R_purify(ispin)%matrix, error=error) 02375 02376 ! *** Construct Matrix M for Hadamard Product 02377 CALL cp_fm_set_all(admm_env%M_purify(ispin)%matrix,0.0_dp,error=error) 02378 pole = 0.0_dp 02379 DO i=1,nao_aux_fit 02380 pole = Heaviside(admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(i)-0.5_dp) 02381 CALL cp_fm_set_element(admm_env%M_purify(ispin)%matrix,i,i,pole,error) 02382 END DO 02383 CALL cp_fm_upper_to_full(admm_env%M_purify(ispin)%matrix,admm_env%work_aux_aux,error=error) 02384 02385 CALL copy_dbcsr_to_fm(density_matrix,admm_env%work_aux_aux3,error) 02386 CALL cp_fm_upper_to_full(admm_env%work_aux_aux3,admm_env%work_aux_aux,error=error) 02387 02388 ! ** S^(-1)*R 02389 CALL cp_fm_gemm('N','N',nao_aux_fit,nao_aux_fit,nao_aux_fit,& 02390 1.0_dp,admm_env%S_inv,admm_env%R_purify(ispin)%matrix,0.0_dp,& 02391 admm_env%work_aux_aux,error) 02392 ! ** S^(-1)*R*M 02393 CALL cp_fm_gemm('N','N',nao_aux_fit,nao_aux_fit,nao_aux_fit,& 02394 1.0_dp,admm_env%work_aux_aux,admm_env%M_purify(ispin)%matrix,0.0_dp,& 02395 admm_env%work_aux_aux2,error) 02396 ! ** S^(-1)*R*M*R^T*S^(-1) 02397 CALL cp_fm_gemm('N','T',nao_aux_fit,nao_aux_fit,nao_aux_fit,& 02398 1.0_dp,admm_env%work_aux_aux2,admm_env%work_aux_aux,0.0_dp,& 02399 admm_env%work_aux_aux3,error) 02400 02401 CALL copy_fm_to_dbcsr(admm_env%work_aux_aux3, density_matrix,keep_sparsity=.TRUE.,error=error) 02402 02403 IF( nspins == 1 ) THEN 02404 CALL cp_dbcsr_scale(density_matrix, 2.0_dp, error=error) 02405 END IF 02406 02407 CALL timestop(handle) 02408 02409 END SUBROUTINE purify_density_matrix_cauchy 02410 02411 02412 FUNCTION delta(x) 02413 REAL(KIND=dp), INTENT(IN) :: x 02414 REAL(KIND=dp) :: delta 02415 02416 IF( x == 0.0_dp) THEN 02417 delta = 1.0_dp 02418 ELSE 02419 delta = 0.0_dp 02420 END IF 02421 02422 END FUNCTION delta 02423 02424 FUNCTION Heaviside(x) 02425 REAL(KIND=dp), INTENT(IN) :: x 02426 REAL(KIND=dp) :: Heaviside 02427 02428 IF( x < 0.0_dp ) THEN 02429 Heaviside = 0.0_dp 02430 ELSE 02431 Heaviside = 1.0_dp 02432 END IF 02433 END FUNCTION Heaviside 02434 02435 END MODULE admm_methods
1.7.3