|
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 ! ***************************************************************************** 00016 MODULE qs_scf_block_davidson 00017 00018 USE array_types, ONLY: array_i1d_obj,& 00019 array_release 00020 USE cp_dbcsr_interface, ONLY: & 00021 cp_create_bl_distribution, cp_dbcsr_add, cp_dbcsr_col_block_sizes, & 00022 cp_dbcsr_copy, cp_dbcsr_create, cp_dbcsr_distribution, & 00023 cp_dbcsr_distribution_release, cp_dbcsr_get_diag, cp_dbcsr_get_info, & 00024 cp_dbcsr_init_p, cp_dbcsr_multiply, cp_dbcsr_norm, cp_dbcsr_release_p, & 00025 cp_dbcsr_row_block_sizes, cp_dbcsr_scale_by_vector 00026 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 00027 copy_fm_to_dbcsr,& 00028 cp_dbcsr_sm_fm_multiply 00029 USE cp_dbcsr_types, ONLY: cp_dbcsr_type 00030 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,& 00031 cp_fm_gemm,& 00032 cp_fm_scale_and_add,& 00033 cp_fm_symm,& 00034 cp_fm_transpose,& 00035 cp_fm_triangular_invert,& 00036 cp_fm_upper_to_full 00037 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,& 00038 cp_fm_cholesky_restore 00039 USE cp_fm_diag, ONLY: choose_eigv_solver,& 00040 cp_fm_power 00041 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 00042 cp_fm_struct_release,& 00043 cp_fm_struct_type 00044 USE cp_fm_types, ONLY: & 00045 cp_fm_create, cp_fm_get_diag, cp_fm_get_info, cp_fm_release, & 00046 cp_fm_set_all, cp_fm_to_fm, cp_fm_to_fm_submat, cp_fm_type, & 00047 cp_fm_vectorsnorm 00048 USE dbcsr_methods, ONLY: dbcsr_distribution_mp,& 00049 dbcsr_distribution_new,& 00050 dbcsr_distribution_row_dist,& 00051 dbcsr_mp_npcols,& 00052 dbcsr_mp_nprows 00053 USE dbcsr_types, ONLY: dbcsr_distribution_obj,& 00054 dbcsr_norm_column,& 00055 dbcsr_type_no_symmetry,& 00056 dbcsr_type_real_default,& 00057 dbcsr_type_symmetric 00058 USE kinds, ONLY: dp 00059 USE machine, ONLY: m_walltime 00060 USE message_passing, ONLY: mp_sum 00061 USE preconditioner, ONLY: apply_preconditioner 00062 USE preconditioner_types, ONLY: preconditioner_type 00063 USE qs_block_davidson_types, ONLY: davidson_type 00064 USE qs_mo_types, ONLY: get_mo_set,& 00065 mo_set_type 00066 USE timings, ONLY: timeset,& 00067 timestop 00068 #include "cp_common_uses.h" 00069 00070 IMPLICIT NONE 00071 PRIVATE 00072 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_block_davidson' 00073 00074 PUBLIC :: generate_extended_space, generate_extended_space_sparse 00075 00076 CONTAINS 00077 00078 ! ***************************************************************************** 00079 00080 SUBROUTINE ritz_coefficients(bdav_env,mo_coeff,matrix_sc,matrix_hc,ritz_coeff,error) 00081 00082 TYPE(davidson_type) :: bdav_env 00083 TYPE(cp_fm_type), POINTER :: mo_coeff, matrix_sc, matrix_hc 00084 REAL(dp), DIMENSION(:) :: ritz_coeff 00085 TYPE(cp_error_type), INTENT(inout) :: error 00086 00087 CHARACTER(len=*), PARAMETER :: routineN = 'ritz_coefficients', 00088 routineP = moduleN//':'//routineN 00089 00090 INTEGER :: handle, i, istat, nao, nmo 00091 LOGICAL :: failure 00092 REAL(dp), ALLOCATABLE, DIMENSION(:) :: chc_diag, csc_diag 00093 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 00094 TYPE(cp_fm_type), POINTER :: block_mat, matrix_tmp 00095 00096 failure=.FALSE. 00097 00098 CALL timeset(routineN,handle) 00099 00100 NULLIFY(block_mat,fm_struct_tmp,matrix_tmp) 00101 CALL cp_fm_get_info(mo_coeff,nrow_global=nao,ncol_global=nmo,error=error) 00102 00103 ALLOCATE(csc_diag(nmo),STAT=istat) 00104 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00105 00106 ALLOCATE(chc_diag(nmo),STAT=istat) 00107 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00108 00109 ! storage matrix of size mos x mos, only the diagonal elements are used 00110 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=nmo,ncol_global=nmo, & 00111 context=mo_coeff%matrix_struct%context, & 00112 para_env=mo_coeff%matrix_struct%para_env,error=error) 00113 CALL cp_fm_create(matrix_tmp,fm_struct_tmp,name="matrix_tmp",error=error) 00114 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00115 00116 ! since we only use diagonal elements this is a bit of a waste 00117 ! compute CSC 00118 ! CALL cp_fm_gemm('T','N',nmo,nmo,nao,1.0_dp,mo_coeff,matrix_sc,0.0_dp,matrix_tmp,error=error) 00119 ! CALL cp_fm_get_diag(matrix_tmp,csc_diag,error=error) 00120 ! set the top left part of S[C,Z] block matrix CSC 00121 block_mat => bdav_env%S_block_mat 00122 CALL cp_fm_set_all(block_mat,0.0_dp, 1.0_dp,error=error) 00123 ! CALL cp_fm_to_fm_submat(matrix_tmp,block_mat,nmo,nmo,1,1,1,1,error=error) 00124 00125 ! compute CHC 00126 CALL cp_fm_gemm('T','N',nmo,nmo,nao,1.0_dp,mo_coeff,matrix_hc,0.0_dp,matrix_tmp,error=error) 00127 CALL cp_fm_get_diag(matrix_tmp,chc_diag,error=error) 00128 ! set the top left part of H[C,Z] block matrix CHC 00129 block_mat => bdav_env%H_block_mat 00130 CALL cp_fm_to_fm_submat(matrix_tmp,block_mat,nmo,nmo,1,1,1,1,error=error) 00131 00132 DO i=1,nmo 00133 ! IF(ABS(csc_diag(i))>EPSILON(0.0_dp)) THEN 00134 ritz_coeff(i) = chc_diag(i)!/csc_diag(i) 00135 ! END IF 00136 END DO 00137 CALL cp_fm_release(matrix_tmp,error=error) 00138 00139 CALL timestop(handle) 00140 00141 END SUBROUTINE ritz_coefficients 00142 ! +++++++++++++++++++++++++++++ 00143 00144 SUBROUTINE generate_extended_space(bdav_env,mo_set,matrix_h,matrix_s,output_unit,& 00145 preconditioner,error) 00146 00147 TYPE(davidson_type) :: bdav_env 00148 TYPE(mo_set_type), POINTER :: mo_set 00149 TYPE(cp_dbcsr_type), POINTER :: matrix_h, matrix_s 00150 INTEGER, INTENT(IN) :: output_unit 00151 TYPE(preconditioner_type), OPTIONAL, 00152 POINTER :: preconditioner 00153 TYPE(cp_error_type), INTENT(inout) :: error 00154 00155 CHARACTER(len=*), PARAMETER :: routineN = 'generate_extended_space', 00156 routineP = moduleN//':'//routineN 00157 00158 INTEGER :: handle, homo, i_first, i_last, imo, istat, iter, j, jj, 00159 max_iter, n, nao, nmat, nmat2, nmo, nmo_converged, nmo_not_converged, 00160 nset, nset_conv, nset_not_conv 00161 INTEGER, ALLOCATABLE, DIMENSION(:) :: iconv, inotconv 00162 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: iconv_set, inotconv_set 00163 LOGICAL :: converged, 00164 do_apply_preconditioner, 00165 failure 00166 REAL(dp) :: lambda, max_norm, min_norm, 00167 t1, t2 00168 REAL(dp), ALLOCATABLE, DIMENSION(:) :: ritz_coeff, vnorm 00169 REAL(dp), DIMENSION(:), POINTER :: eig_not_conv, eigenvalues, 00170 evals 00171 TYPE(cp_dbcsr_type), POINTER :: mo_coeff_b 00172 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 00173 TYPE(cp_fm_type), POINTER :: c_conv, c_notconv, c_out, c_pz, c_z, 00174 h_block, h_fm, m_hc, m_sc, m_tmp, mo_coeff, mt_tmp, s_block, s_fm, 00175 v_block, w_block 00176 00177 failure=.FALSE. 00178 00179 CALL timeset(routineN,handle) 00180 00181 NULLIFY(mo_coeff,mo_coeff_b,eigenvalues) 00182 00183 do_apply_preconditioner = .FALSE. 00184 IF(PRESENT(preconditioner)) do_apply_preconditioner=.TRUE. 00185 CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff,mo_coeff_b=mo_coeff_b,eigenvalues=eigenvalues,& 00186 nao=nao,nmo=nmo, homo=homo) 00187 IF(do_apply_preconditioner) THEN 00188 max_iter = bdav_env%max_iter 00189 ELSE 00190 max_iter = 1 00191 END IF 00192 00193 NULLIFY(c_conv, c_notconv, c_out, c_z, c_pz, m_hc, m_sc, m_tmp, mt_tmp) 00194 NULLIFY(h_block, h_fm, s_block, s_fm, v_block, w_block) 00195 NULLIFY(evals,eig_not_conv ) 00196 t1 = m_walltime() 00197 IF (output_unit > 0) THEN 00198 WRITE(output_unit,"(T15,A,T23,A,T36,A,T49,A,T60,A,/,T8,A)") & 00199 " Cycle ", " conv. MOS ", " B2MAX ", " B2MIN ", " Time", REPEAT("-",60) 00200 END IF 00201 00202 ALLOCATE(iconv(nmo), STAT=istat) 00203 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00204 ALLOCATE(inotconv(nmo), STAT=istat) 00205 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00206 ALLOCATE(ritz_coeff(nmo),STAT=istat) 00207 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00208 ALLOCATE(vnorm(nmo),STAT=istat) 00209 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00210 00211 converged=.FALSE. 00212 DO iter = 1, max_iter 00213 00214 ! compute Ritz values 00215 ritz_coeff=0.0_dp 00216 CALL cp_fm_create(m_hc,mo_coeff%matrix_struct,name="hc",error=error) 00217 CALL cp_dbcsr_sm_fm_multiply(matrix_h,mo_coeff,m_hc,nmo,error=error) 00218 CALL cp_fm_create(m_sc,mo_coeff%matrix_struct,name="sc",error=error) 00219 CALL cp_dbcsr_sm_fm_multiply(matrix_s,mo_coeff,m_sc,nmo,error=error) 00220 00221 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=nmo,ncol_global=nmo, & 00222 context=mo_coeff%matrix_struct%context, & 00223 para_env=mo_coeff%matrix_struct%para_env,error=error) 00224 CALL cp_fm_create(m_tmp,fm_struct_tmp,name="matrix_tmp",error=error) 00225 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00226 00227 CALL cp_fm_gemm('T','N',nmo,nmo,nao,1.0_dp,mo_coeff,m_hc,0.0_dp,m_tmp,error=error) 00228 CALL cp_fm_get_diag(m_tmp,ritz_coeff,error=error) 00229 CALL cp_fm_release(m_tmp,error=error) 00230 00231 ! Check for converged eigenvectors 00232 ! CALL cp_fm_create(c_z,mo_coeff%matrix_struct,name="tmp",error=error) 00233 c_z => bdav_env%matrix_z 00234 c_pz => bdav_env%matrix_pz 00235 CALL cp_fm_to_fm(m_sc,c_z,error=error) 00236 CALL cp_fm_column_scale(c_z,ritz_coeff) 00237 CALL cp_fm_scale_and_add(-1.0_dp,c_z,1.0_dp,m_hc,error=error) 00238 CALL cp_fm_vectorsnorm(c_z,vnorm,error=error) 00239 00240 nmo_converged = 0 00241 nmo_not_converged = 0 00242 max_norm = 0.0_dp 00243 min_norm = 1.e10_dp 00244 DO imo = 1,nmo 00245 max_norm = MAX(max_norm,vnorm(imo)) 00246 min_norm = MIN(min_norm,vnorm(imo)) 00247 END DO 00248 iconv = 0 00249 inotconv = 0 00250 DO imo = 1,nmo 00251 IF(vnorm(imo) <= bdav_env%eps_iter ) THEN 00252 nmo_converged = nmo_converged + 1 00253 iconv(nmo_converged)=imo 00254 ELSE 00255 nmo_not_converged = nmo_not_converged + 1 00256 inotconv(nmo_not_converged)=imo 00257 END IF 00258 END DO 00259 00260 IF(nmo_converged>0) THEN 00261 ALLOCATE(iconv_set(nmo_converged,2), STAT=istat) 00262 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00263 ALLOCATE(inotconv_set(nmo_not_converged,2), STAT=istat) 00264 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00265 i_last = iconv(1) 00266 nset = 0 00267 DO j = 1,nmo_converged 00268 imo=iconv(j) 00269 00270 IF(imo==i_last+1) THEN 00271 i_last=imo 00272 iconv_set(nset,2)=imo 00273 ELSE 00274 i_last=imo 00275 nset=nset+1 00276 iconv_set(nset,1)=imo 00277 iconv_set(nset,2)=imo 00278 END IF 00279 END DO 00280 nset_conv = nset 00281 00282 i_last = inotconv(1) 00283 nset = 0 00284 DO j = 1,nmo_not_converged 00285 imo=inotconv(j) 00286 00287 IF(imo==i_last+1) THEN 00288 i_last=imo 00289 inotconv_set(nset,2)=imo 00290 ELSE 00291 i_last=imo 00292 nset=nset+1 00293 inotconv_set(nset,1)=imo 00294 inotconv_set(nset,2)=imo 00295 END IF 00296 END DO 00297 nset_not_conv = nset 00298 CALL cp_fm_release(m_sc,error=error) 00299 CALL cp_fm_release(m_hc,error=error) 00300 NULLIFY(c_z, c_pz) 00301 END IF 00302 00303 00304 IF(REAL(nmo_converged,dp)/REAL(nmo,dp)>bdav_env%conv_percent) THEN 00305 converged=.TRUE. 00306 DEALLOCATE(iconv_set,STAT=istat) 00307 DEALLOCATE (inotconv_set,STAT=istat) 00308 t2 = m_walltime() 00309 IF (output_unit > 0) THEN 00310 WRITE(output_unit,'(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)')& 00311 iter, nmo_converged, max_norm, min_norm, t2-t1 00312 00313 WRITE(output_unit,*) " Reached convergence in ", iter, & 00314 " Davidson iterations" 00315 END IF 00316 00317 EXIT 00318 END IF 00319 00320 IF(nmo_converged>0) THEN 00321 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=nao,ncol_global=nao, & 00322 context=mo_coeff%matrix_struct%context, & 00323 para_env=mo_coeff%matrix_struct%para_env,error=error) 00324 !allocate h_fm 00325 CALL cp_fm_create(h_fm,fm_struct_tmp,name="matrix_tmp",error=error) 00326 !allocate s_fm 00327 CALL cp_fm_create(s_fm,fm_struct_tmp,name="matrix_tmp",error=error) 00328 !copy matrix_h in h_fm 00329 CALL copy_dbcsr_to_fm(matrix_h,h_fm,error=error) 00330 CALL cp_fm_upper_to_full(h_fm,s_fm,error=error) 00331 00332 !copy matrix_s in s_fm 00333 ! CALL cp_fm_set_all(s_fm,0.0_dp,error=error) 00334 CALL copy_dbcsr_to_fm(matrix_s,s_fm,error=error) 00335 00336 !allocate c_out 00337 CALL cp_fm_create(c_out,fm_struct_tmp,name="matrix_tmp",error=error) 00338 ! set c_out to zero 00339 CALL cp_fm_set_all(c_out,0.0_dp,error=error) 00340 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00341 00342 !allocate c_conv 00343 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=nao,ncol_global=nmo_converged, & 00344 context=mo_coeff%matrix_struct%context, & 00345 para_env=mo_coeff%matrix_struct%para_env,error=error) 00346 CALL cp_fm_create(c_conv,fm_struct_tmp,name="c_conv",error=error) 00347 CALL cp_fm_set_all(c_conv,0.0_dp,error=error) 00348 !allocate m_tmp 00349 CALL cp_fm_create(m_tmp,fm_struct_tmp,name="m_tmp_nxmc",error=error) 00350 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00351 END IF 00352 00353 !allocate c_notconv 00354 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=nao,ncol_global=nmo_not_converged, & 00355 context=mo_coeff%matrix_struct%context, & 00356 para_env=mo_coeff%matrix_struct%para_env,error=error) 00357 CALL cp_fm_create(c_notconv,fm_struct_tmp,name="c_notconv",error=error) 00358 CALL cp_fm_set_all(c_notconv,0.0_dp,error=error) 00359 IF(nmo_converged>0) THEN 00360 CALL cp_fm_create(m_hc,fm_struct_tmp,name="m_hc",error=error) 00361 CALL cp_fm_create(m_sc,fm_struct_tmp,name="m_sc",error=error) 00362 !allocate c_z 00363 CALL cp_fm_create(c_z,fm_struct_tmp,name="c_z",error=error) 00364 CALL cp_fm_create(c_pz,fm_struct_tmp,name="c_pz",error=error) 00365 CALL cp_fm_set_all(c_z,0.0_dp,error=error) 00366 00367 ! sum contributions to c_out 00368 jj=1 00369 DO j=1,nset_conv 00370 i_first=iconv_set(j,1) 00371 i_last=iconv_set(j,2) 00372 n=i_last-i_first+1 00373 CALL cp_fm_to_fm_submat(mo_coeff,c_conv,nao,n,1,i_first,1,jj,error=error) 00374 jj=jj+n 00375 END DO 00376 CALL cp_fm_symm('L','U',nao,nmo_converged,1.0_dp,s_fm,c_conv,0.0_dp,m_tmp,error=error) 00377 CALL cp_fm_gemm('N','T',nao,nao,nmo_converged,1.0_dp,m_tmp,m_tmp,0.0_dp,c_out,error=error) 00378 00379 ! project c_out out of H 00380 lambda = 100.0_dp*ABS(eigenvalues(homo)) 00381 CALL cp_fm_scale_and_add(lambda,c_out,1.0_dp,h_fm,error=error) 00382 CALL cp_fm_release(m_tmp,error=error) 00383 CALL cp_fm_release(h_fm,error=error) 00384 00385 END IF 00386 00387 !allocate m_tmp 00388 CALL cp_fm_create(m_tmp,fm_struct_tmp,name="m_tmp_nxm",error=error) 00389 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00390 IF(nmo_converged>0) THEN 00391 ALLOCATE(eig_not_conv(nmo_not_converged)) 00392 jj=1 00393 DO j=1,nset_not_conv 00394 i_first=inotconv_set(j,1) 00395 i_last=inotconv_set(j,2) 00396 n=i_last-i_first+1 00397 CALL cp_fm_to_fm_submat(mo_coeff,c_notconv,nao,n,1,i_first,1,jj,error=error) 00398 eig_not_conv(jj:jj+n-1) = ritz_coeff(i_first:i_last) 00399 jj=jj+n 00400 END DO 00401 CALL cp_fm_gemm('N','N',nao,nmo_not_converged,nao,1.0_dp,c_out,c_notconv,0.0_dp,m_hc,error=error) 00402 CALL cp_fm_symm('L','U',nao,nmo_not_converged,1.0_dp,s_fm,c_notconv,0.0_dp,m_sc,error=error) 00403 ! extend suspace using only the not converged vectors 00404 CALL cp_fm_to_fm(m_sc,m_tmp,error=error) 00405 CALL cp_fm_column_scale(m_tmp,eig_not_conv) 00406 CALL cp_fm_scale_and_add(-1.0_dp,m_tmp,1.0_dp,m_hc,error=error) 00407 DEALLOCATE(eig_not_conv) 00408 CALL cp_fm_to_fm(m_tmp,c_z,error=error) 00409 ELSE 00410 CALL cp_fm_to_fm(mo_coeff,c_notconv,error=error) 00411 END IF 00412 00413 !preconditioner 00414 IF(do_apply_preconditioner) THEN 00415 IF (preconditioner%in_use/=0) THEN 00416 CALL apply_preconditioner(preconditioner,c_z,c_pz,error=error) 00417 ELSE 00418 CALL cp_fm_to_fm(c_z,c_pz,error=error) 00419 ENDIF 00420 ELSE 00421 CALL cp_fm_to_fm(c_z,c_pz,error=error) 00422 END IF 00423 CALL cp_fm_release(m_tmp,error=error) 00424 00425 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=nmo_not_converged,ncol_global=nmo_not_converged, & 00426 context=mo_coeff%matrix_struct%context, & 00427 para_env=mo_coeff%matrix_struct%para_env,error=error) 00428 00429 CALL cp_fm_create(m_tmp,fm_struct_tmp,name="m_tmp_mxm",error=error) 00430 CALL cp_fm_create(mt_tmp,fm_struct_tmp,name="mt_tmp_mxm",error=error) 00431 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00432 00433 nmat = nmo_not_converged 00434 nmat2 = 2* nmo_not_converged 00435 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=nmat2,ncol_global=nmat2, & 00436 context=mo_coeff%matrix_struct%context, & 00437 para_env=mo_coeff%matrix_struct%para_env,error=error) 00438 00439 CALL cp_fm_create(s_block,fm_struct_tmp,name="sb",error=error) 00440 CALL cp_fm_create(h_block,fm_struct_tmp,name="hb",error=error) 00441 CALL cp_fm_create(v_block,fm_struct_tmp,name="vb",error=error) 00442 CALL cp_fm_create(w_block,fm_struct_tmp,name="wb",error=error) 00443 ALLOCATE(evals(nmat2)) 00444 00445 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00446 00447 ! compute CSC 00448 CALL cp_fm_set_all(s_block,0.0_dp, 1.0_dp,error=error) 00449 00450 ! compute CHC 00451 CALL cp_fm_gemm('T','N',nmat,nmat,nao,1.0_dp,c_notconv,m_hc,0.0_dp,m_tmp,error=error) 00452 CALL cp_fm_to_fm_submat(m_tmp,h_block,nmat,nmat,1,1,1,1,error=error) 00453 00454 ! compute ZSC 00455 CALL cp_fm_gemm('T','N',nmat,nmat,nao,1.0_dp,c_pz,m_sc,0.0_dp,m_tmp,error=error) 00456 CALL cp_fm_to_fm_submat(m_tmp,s_block,nmat,nmat,1,1,1+nmat,1,error=error) 00457 CALL cp_fm_transpose(m_tmp,mt_tmp,error=error) 00458 CALL cp_fm_to_fm_submat(mt_tmp,s_block,nmat,nmat,1,1,1,1+nmat,error=error) 00459 ! compute ZHC 00460 CALL cp_fm_gemm('T','N',nmat,nmat,nao,1.0_dp,c_pz,m_hc,0.0_dp,m_tmp,error=error) 00461 CALL cp_fm_to_fm_submat(m_tmp,h_block,nmat,nmat,1,1,1+nmat,1,error=error) 00462 CALL cp_fm_transpose(m_tmp,mt_tmp,error=error) 00463 CALL cp_fm_to_fm_submat(mt_tmp,h_block,nmat,nmat,1,1,1,1+nmat,error=error) 00464 00465 CALL cp_fm_release(mt_tmp,error=error) 00466 00467 ! reuse m_sc and m_hc to computr HZ and SZ 00468 IF(nmo_converged>0) THEN 00469 CALL cp_fm_gemm('N','N',nao,nmat,nao,1.0_dp,c_out,c_pz,0.0_dp,m_hc,error=error) 00470 CALL cp_fm_symm('L','U',nao,nmo_not_converged,1.0_dp,s_fm,c_pz,0.0_dp,m_sc,error=error) 00471 00472 CALL cp_fm_release(c_out,error=error) 00473 CALL cp_fm_release(c_conv,error=error) 00474 CALL cp_fm_release(s_fm,error=error) 00475 ELSE 00476 CALL cp_dbcsr_sm_fm_multiply(matrix_h,c_pz,m_hc,nmo,error=error) 00477 CALL cp_dbcsr_sm_fm_multiply(matrix_s,c_pz,m_sc,nmo,error=error) 00478 END IF 00479 00480 ! compute ZSZ 00481 CALL cp_fm_gemm('T','N',nmat,nmat,nao,1.0_dp,c_pz,m_sc,0.0_dp,m_tmp,error=error) 00482 CALL cp_fm_to_fm_submat(m_tmp,s_block,nmat,nmat,1,1,1+nmat,1+nmat,error=error) 00483 ! compute ZHZ 00484 CALL cp_fm_gemm('T','N',nmat,nmat,nao,1.0_dp,c_pz,m_hc,0.0_dp,m_tmp,error=error) 00485 CALL cp_fm_to_fm_submat(m_tmp,h_block,nmat,nmat,1,1,1+nmat,1+nmat,error=error) 00486 00487 CALL cp_fm_release(m_sc,error=error) 00488 00489 ! solution of the reduced eigenvalues problem 00490 CALL reduce_extended_space(s_block,h_block,v_block,w_block, evals, nmat2, error=error) 00491 00492 ! extract egenvectors 00493 CALL cp_fm_to_fm_submat(v_block,m_tmp, nmat,nmat,1,1,1,1,error=error) 00494 CALL cp_fm_gemm('N','N',nao,nmat,nmat,1.0_dp,c_notconv,m_tmp,0.0_dp,m_hc,error=error) 00495 CALL cp_fm_to_fm_submat(v_block,m_tmp, nmat,nmat,1+nmat,1,1,1,error=error) 00496 CALL cp_fm_gemm('N','N',nao,nmat,nmat,1.0_dp,c_pz,m_tmp,1.0_dp,m_hc,error=error) 00497 00498 CALL cp_fm_release(m_tmp,error=error) 00499 00500 CALL cp_fm_release(c_notconv,error=error) 00501 CALL cp_fm_release(s_block,error=error) 00502 CALL cp_fm_release(h_block,error=error) 00503 CALL cp_fm_release(w_block,error=error) 00504 CALL cp_fm_release(v_block,error=error) 00505 00506 IF(nmo_converged>0) THEN 00507 CALL cp_fm_release(c_z,error=error) 00508 CALL cp_fm_release(c_pz,error=error) 00509 jj=1 00510 DO j=1,nset_not_conv 00511 i_first=inotconv_set(j,1) 00512 i_last=inotconv_set(j,2) 00513 n=i_last-i_first+1 00514 CALL cp_fm_to_fm_submat(m_hc,mo_coeff,nao,n,1,jj,1,i_first,error=error) 00515 eigenvalues(i_first:i_last) = evals(jj:jj+n-1) 00516 jj=jj+n 00517 END DO 00518 DEALLOCATE(iconv_set,STAT=istat) 00519 DEALLOCATE (inotconv_set,STAT=istat) 00520 ELSE 00521 CALL cp_fm_to_fm(m_hc,mo_coeff,error=error) 00522 eigenvalues(1:nmo) = evals(1:nmo) 00523 END IF 00524 DEALLOCATE(evals,STAT=istat) 00525 00526 CALL cp_fm_release(m_hc,error=error) 00527 00528 CALL copy_fm_to_dbcsr(mo_coeff,mo_coeff_b,error=error)!fm->dbcsr 00529 00530 t2 = m_walltime() 00531 IF (output_unit > 0) THEN 00532 WRITE(output_unit,'(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)')& 00533 iter, nmo_converged, max_norm, min_norm, t2-t1 00534 END IF 00535 t1=m_walltime() 00536 00537 END DO ! iter 00538 00539 DEALLOCATE(iconv,STAT=istat) 00540 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00541 DEALLOCATE(inotconv,STAT=istat) 00542 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00543 DEALLOCATE(ritz_coeff,STAT=istat) 00544 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00545 DEALLOCATE(vnorm,STAT=istat) 00546 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00547 00548 CALL timestop(handle) 00549 END SUBROUTINE generate_extended_space 00550 00551 ! ***************************************************************************** 00552 SUBROUTINE generate_extended_space_sparse(bdav_env,mo_set,matrix_h,matrix_s,output_unit,& 00553 preconditioner,error) 00554 00555 TYPE(davidson_type) :: bdav_env 00556 TYPE(mo_set_type), POINTER :: mo_set 00557 TYPE(cp_dbcsr_type), POINTER :: matrix_h, matrix_s 00558 INTEGER, INTENT(IN) :: output_unit 00559 TYPE(preconditioner_type), OPTIONAL, 00560 POINTER :: preconditioner 00561 TYPE(cp_error_type), INTENT(inout) :: error 00562 00563 CHARACTER(len=*), PARAMETER :: 00564 routineN = 'generate_extended_space_sparse', 00565 routineP = moduleN//':'//routineN 00566 00567 INTEGER :: handle, homo, i_first, i_last, imo, istat, iter, j, jj, k, 00568 max_iter, n, nao, nmat, nmat2, nmo, nmo_converged, nmo_not_converged, 00569 nset, nset_conv, nset_not_conv 00570 INTEGER, ALLOCATABLE, DIMENSION(:) :: iconv, inotconv 00571 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: iconv_set, inotconv_set 00572 LOGICAL :: converged, 00573 do_apply_preconditioner, 00574 failure 00575 REAL(dp) :: lambda, max_norm, min_norm, 00576 t1, t2 00577 REAL(dp), ALLOCATABLE, DIMENSION(:) :: eig_not_conv, evals, 00578 ritz_coeff, vnorm 00579 REAL(dp), DIMENSION(:), POINTER :: eigenvalues 00580 TYPE(array_i1d_obj) :: col_blk_size, col_dist, 00581 row_blk_size, row_dist 00582 TYPE(cp_dbcsr_type), POINTER :: c_out, matrix_hc, matrix_mm, matrix_pz, 00583 matrix_sc, matrix_z, mo_coeff_b, mo_conv, mo_notconv, smo_conv 00584 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 00585 TYPE(cp_fm_type), POINTER :: h_block, matrix_mm_fm, matrix_mmt_fm, 00586 matrix_nm_fm, matrix_z_fm, mo_coeff, mo_conv_fm, mo_notconv_fm, 00587 s_block, v_block, w_block 00588 TYPE(dbcsr_distribution_obj) :: dist 00589 00590 failure=.FALSE. 00591 CALL timeset(routineN,handle) 00592 00593 do_apply_preconditioner = .FALSE. 00594 IF(PRESENT(preconditioner)) do_apply_preconditioner=.TRUE. 00595 00596 NULLIFY(mo_coeff,mo_coeff_b, matrix_hc, matrix_sc, matrix_z, matrix_pz, matrix_mm) 00597 NULLIFY(mo_conv_fm, mo_notconv_fm, mo_conv, mo_notconv, smo_conv,c_out) 00598 NULLIFY(matrix_mm_fm, matrix_mmt_fm, mo_coeff, matrix_nm_fm, matrix_z_fm) 00599 NULLIFY(h_block,s_block,v_block,w_block) 00600 NULLIFY(fm_struct_tmp) 00601 CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b,& 00602 eigenvalues=eigenvalues,homo=homo,nao=nao,nmo=nmo) 00603 IF(do_apply_preconditioner) THEN 00604 max_iter = bdav_env%max_iter 00605 ELSE 00606 max_iter = 1 00607 END IF 00608 00609 t1 = m_walltime() 00610 IF (output_unit > 0) THEN 00611 WRITE(output_unit,"(T15,A,T23,A,T36,A,T49,A,T60,A,/,T8,A)") & 00612 " Cycle ", " conv. MOS ", " B2MAX ", " B2MIN ", " Time", REPEAT("-",60) 00613 END IF 00614 00615 00616 ! Allocate array for Ritz values 00617 ALLOCATE(ritz_coeff(nmo),STAT=istat) 00618 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00619 ALLOCATE(iconv(nmo), STAT=istat) 00620 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00621 ALLOCATE(inotconv(nmo), STAT=istat) 00622 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00623 ALLOCATE(vnorm(nmo),STAT=istat) 00624 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00625 00626 converged=.FALSE. 00627 DO iter = 1, max_iter 00628 NULLIFY(c_out, mo_conv, mo_notconv_fm, mo_notconv) 00629 ! Prepare HC and SC, using mo_coeff_b (sparse), these are still sparse 00630 CALL cp_dbcsr_init_p(matrix_hc,error=error) 00631 CALL cp_dbcsr_create(matrix_hc,"matrix_hc",cp_dbcsr_distribution(mo_coeff_b),& 00632 dbcsr_type_no_symmetry,cp_dbcsr_row_block_sizes(mo_coeff_b),& 00633 cp_dbcsr_col_block_sizes(mo_coeff_b),0,0,dbcsr_type_real_default,error=error) 00634 CALL cp_dbcsr_init_p(matrix_sc,error=error) 00635 CALL cp_dbcsr_create(matrix_sc,"matrix_sc",cp_dbcsr_distribution(mo_coeff_b),& 00636 dbcsr_type_no_symmetry,cp_dbcsr_row_block_sizes(mo_coeff_b),& 00637 cp_dbcsr_col_block_sizes(mo_coeff_b),0,0,dbcsr_type_real_default,error=error) 00638 00639 CALL cp_dbcsr_get_info(mo_coeff_b,nfullrows_total=n,nfullcols_total=k) 00640 CALL cp_dbcsr_multiply('n','n',1.0_dp,matrix_h,mo_coeff_b,0.0_dp,matrix_hc,last_column=k,error=error) 00641 CALL cp_dbcsr_multiply('n','n',1.0_dp,matrix_s,mo_coeff_b,0.0_dp,matrix_sc,last_column=k,error=error) 00642 00643 ! compute Ritz values 00644 ritz_coeff=0.0_dp 00645 ! Allocate Sparse matrices: nmoxnmo 00646 ! matrix_mm 00647 CALL cp_create_bl_distribution (col_dist, col_blk_size, nmo, & 00648 dbcsr_mp_npcols(dbcsr_distribution_mp(cp_dbcsr_distribution(mo_coeff_b)))) 00649 CALL cp_create_bl_distribution (row_dist, row_blk_size, nmo, & 00650 dbcsr_mp_nprows(dbcsr_distribution_mp(cp_dbcsr_distribution(mo_coeff_b)))) 00651 CALL dbcsr_distribution_new (dist, dbcsr_distribution_mp(cp_dbcsr_distribution(mo_coeff_b)),& 00652 row_dist,col_dist) 00653 CALL cp_dbcsr_init_p(matrix_mm,error=error) 00654 CALL cp_dbcsr_create(matrix_mm,"matrix_mm",dist,dbcsr_type_no_symmetry,& 00655 row_blk_size,col_blk_size,0,0,dbcsr_type_real_default,error=error) 00656 CALL cp_dbcsr_distribution_release (dist) 00657 CALL array_release (col_blk_size) 00658 CALL array_release (col_dist) 00659 CALL array_release (row_blk_size) 00660 CALL array_release (row_dist) 00661 00662 00663 CALL cp_dbcsr_multiply('t','n',1.0_dp,mo_coeff_b,matrix_hc,0.0_dp,matrix_mm, last_column=k, error=error) 00664 CALL cp_dbcsr_get_diag(matrix_mm,ritz_coeff,error=error) 00665 CALL mp_sum(ritz_coeff,mo_coeff%matrix_struct%para_env%group) 00666 00667 00668 ! extended subspace P Z = P [H - theta S]C this ia another matrix of type and size as mo_coeff_b 00669 CALL cp_dbcsr_init_p(matrix_z,error=error) 00670 CALL cp_dbcsr_create(matrix_z,"matrix_z",cp_dbcsr_distribution(mo_coeff_b),& 00671 dbcsr_type_no_symmetry,cp_dbcsr_row_block_sizes(mo_coeff_b),& 00672 cp_dbcsr_col_block_sizes(mo_coeff_b),0,0,dbcsr_type_real_default,error=error) 00673 CALL cp_dbcsr_copy(matrix_z,matrix_sc,error=error) 00674 CALL cp_dbcsr_scale_by_vector(matrix_z,ritz_coeff,side='right',error=error) 00675 CALL cp_dbcsr_add(matrix_z,matrix_hc,-1.0_dp,1.0_dp,error=error) 00676 00677 ! Check for converged eigenvectors 00678 vnorm =0.0_dp 00679 CALL cp_dbcsr_norm(matrix_z,which_norm=dbcsr_norm_column,norm_vector=vnorm,error=error) 00680 nmo_converged = 0 00681 nmo_not_converged = 0 00682 max_norm = 0.0_dp 00683 min_norm = 1.e10_dp 00684 DO imo = 1,nmo 00685 max_norm = MAX(max_norm,vnorm(imo)) 00686 min_norm = MIN(min_norm,vnorm(imo)) 00687 END DO 00688 iconv = 0 00689 inotconv = 0 00690 00691 DO imo = 1,nmo 00692 IF(vnorm(imo) <= bdav_env%eps_iter ) THEN 00693 nmo_converged = nmo_converged + 1 00694 iconv(nmo_converged)=imo 00695 ELSE 00696 nmo_not_converged = nmo_not_converged + 1 00697 inotconv(nmo_not_converged)=imo 00698 END IF 00699 END DO 00700 00701 00702 IF(nmo_converged>0) THEN 00703 ALLOCATE(iconv_set(nmo_converged,2), STAT=istat) 00704 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00705 ALLOCATE(inotconv_set(nmo_not_converged,2), STAT=istat) 00706 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00707 i_last = iconv(1) 00708 nset = 0 00709 DO j = 1,nmo_converged 00710 imo=iconv(j) 00711 00712 IF(imo==i_last+1) THEN 00713 i_last=imo 00714 iconv_set(nset,2)=imo 00715 ELSE 00716 i_last=imo 00717 nset=nset+1 00718 iconv_set(nset,1)=imo 00719 iconv_set(nset,2)=imo 00720 END IF 00721 END DO 00722 nset_conv = nset 00723 00724 i_last = inotconv(1) 00725 nset = 0 00726 DO j = 1,nmo_not_converged 00727 imo=inotconv(j) 00728 00729 IF(imo==i_last+1) THEN 00730 i_last=imo 00731 inotconv_set(nset,2)=imo 00732 ELSE 00733 i_last=imo 00734 nset=nset+1 00735 inotconv_set(nset,1)=imo 00736 inotconv_set(nset,2)=imo 00737 END IF 00738 END DO 00739 nset_not_conv = nset 00740 00741 CALL cp_dbcsr_release_p(matrix_hc, error=error) 00742 CALL cp_dbcsr_release_p(matrix_sc, error=error) 00743 CALL cp_dbcsr_release_p(matrix_z, error=error) 00744 CALL cp_dbcsr_release_p(matrix_mm, error=error) 00745 END IF 00746 00747 IF(REAL(nmo_converged,dp)/REAL(nmo,dp)>bdav_env%conv_percent) THEN 00748 DEALLOCATE(iconv_set,STAT=istat) 00749 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00750 00751 DEALLOCATE (inotconv_set,STAT=istat) 00752 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00753 00754 converged=.TRUE. 00755 t2 = m_walltime() 00756 IF (output_unit > 0) THEN 00757 WRITE(output_unit,'(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)')& 00758 iter, nmo_converged, max_norm, min_norm, t2-t1 00759 00760 WRITE(output_unit,*) " Reached convergence in ", iter, & 00761 " Davidson iterations" 00762 END IF 00763 00764 EXIT 00765 END IF 00766 00767 IF(nmo_converged>0) THEN 00768 00769 !allocate mo_conv_fm 00770 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=nao,ncol_global=nmo_converged, & 00771 context=mo_coeff%matrix_struct%context, & 00772 para_env=mo_coeff%matrix_struct%para_env,error=error) 00773 CALL cp_fm_create(mo_conv_fm,fm_struct_tmp,name="mo_conv_fm",error=error) 00774 00775 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00776 00777 ! extract mo_conv from mo_coeff full matrix 00778 jj=1 00779 DO j=1,nset_conv 00780 i_first=iconv_set(j,1) 00781 i_last=iconv_set(j,2) 00782 n=i_last-i_first+1 00783 CALL cp_fm_to_fm_submat(mo_coeff,mo_conv_fm,nao,n,1,i_first,1,jj,error=error) 00784 jj=jj+n 00785 END DO 00786 00787 ! allocate c_out sparse matrix, to project out the converged MOS 00788 CALL cp_dbcsr_init_p(c_out,error=error) 00789 CALL cp_dbcsr_create(c_out,"c_out",cp_dbcsr_distribution(matrix_s),& 00790 dbcsr_type_symmetric,cp_dbcsr_row_block_sizes(matrix_s),& 00791 cp_dbcsr_col_block_sizes(matrix_s),0,0,dbcsr_type_real_default,error=error) 00792 00793 ! allocate mo_conv sparse 00794 CALL cp_create_bl_distribution (col_dist, col_blk_size, nmo_converged, & 00795 dbcsr_mp_npcols(dbcsr_distribution_mp(cp_dbcsr_distribution(matrix_s)))) 00796 CALL dbcsr_distribution_new (dist, dbcsr_distribution_mp(cp_dbcsr_distribution(matrix_s)),& 00797 dbcsr_distribution_row_dist(cp_dbcsr_distribution(matrix_s)),col_dist) 00798 CALL cp_dbcsr_init_p(mo_conv,error=error) 00799 CALL cp_dbcsr_create(mo_conv,"mo_conv",dist,dbcsr_type_no_symmetry,& 00800 cp_dbcsr_row_block_sizes(matrix_s),col_blk_size,0,0,dbcsr_type_real_default,error=error) 00801 CALL cp_dbcsr_init_p(smo_conv,error=error) 00802 CALL cp_dbcsr_create(smo_conv,"smo_conv",dist,dbcsr_type_no_symmetry,& 00803 cp_dbcsr_row_block_sizes(matrix_s),col_blk_size,0,0,dbcsr_type_real_default,error=error) 00804 00805 CALL copy_fm_to_dbcsr(mo_conv_fm,mo_conv,error=error)!fm->dbcsr 00806 CALL cp_dbcsr_distribution_release (dist) 00807 CALL array_release (col_blk_size) 00808 CALL array_release (col_dist) 00809 00810 CALL cp_dbcsr_multiply('n','n',1.0_dp,matrix_s,mo_conv,0.0_dp,smo_conv,last_column=nmo_converged,error=error) 00811 CALL cp_dbcsr_multiply('n','t',1.0_dp,smo_conv,smo_conv,0.0_dp,c_out,last_column=nao,error=error) 00812 ! project c_out out of H 00813 lambda = 100.0_dp*ABS(eigenvalues(homo)) 00814 CALL cp_dbcsr_add(c_out,matrix_h,lambda,1.0_dp,error=error) 00815 00816 CALL cp_dbcsr_release_p(mo_conv, error=error) 00817 CALL cp_dbcsr_release_p(smo_conv, error=error) 00818 CALL cp_fm_release(mo_conv_fm, error=error) 00819 00820 00821 !allocate c_notconv_fm 00822 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=nao,ncol_global=nmo_not_converged, & 00823 context=mo_coeff%matrix_struct%context, & 00824 para_env=mo_coeff%matrix_struct%para_env,error=error) 00825 CALL cp_fm_create(mo_notconv_fm,fm_struct_tmp,name="mo_notconv_fm",error=error) 00826 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00827 00828 ! extract mo_notconv from mo_coeff full matrix 00829 ALLOCATE(eig_not_conv(nmo_not_converged),STAT=istat) 00830 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00831 00832 jj=1 00833 DO j=1,nset_not_conv 00834 i_first=inotconv_set(j,1) 00835 i_last=inotconv_set(j,2) 00836 n=i_last-i_first+1 00837 CALL cp_fm_to_fm_submat(mo_coeff,mo_notconv_fm,nao,n,1,i_first,1,jj,error=error) 00838 eig_not_conv(jj:jj+n-1) = ritz_coeff(i_first:i_last) 00839 jj=jj+n 00840 END DO 00841 00842 ! allocate mo_conv sparse 00843 CALL cp_create_bl_distribution (col_dist, col_blk_size, nmo_not_converged, & 00844 dbcsr_mp_npcols(dbcsr_distribution_mp(cp_dbcsr_distribution(matrix_s)))) 00845 CALL dbcsr_distribution_new (dist, dbcsr_distribution_mp(cp_dbcsr_distribution(matrix_s)),& 00846 dbcsr_distribution_row_dist(cp_dbcsr_distribution(matrix_s)),col_dist) 00847 CALL cp_dbcsr_init_p(mo_notconv,error=error) 00848 CALL cp_dbcsr_create(mo_notconv,"mo_notconv",dist,dbcsr_type_no_symmetry,& 00849 cp_dbcsr_row_block_sizes(matrix_s),col_blk_size,0,0,dbcsr_type_real_default,error=error) 00850 00851 CALL cp_dbcsr_init_p(matrix_hc,error=error) 00852 CALL cp_dbcsr_create(matrix_hc,"matrix_hc",dist,dbcsr_type_no_symmetry,& 00853 cp_dbcsr_row_block_sizes(matrix_s),col_blk_size,0,0,dbcsr_type_real_default,error=error) 00854 00855 CALL cp_dbcsr_init_p(matrix_sc,error=error) 00856 CALL cp_dbcsr_create(matrix_sc,"matrix_sc",dist,dbcsr_type_no_symmetry,& 00857 cp_dbcsr_row_block_sizes(matrix_s),col_blk_size,0,0,dbcsr_type_real_default,error=error) 00858 00859 CALL cp_dbcsr_init_p(matrix_z,error=error) 00860 CALL cp_dbcsr_create(matrix_z,"matrix_z",dist,dbcsr_type_no_symmetry,& 00861 cp_dbcsr_row_block_sizes(matrix_s),col_blk_size,0,0,dbcsr_type_real_default,error=error) 00862 00863 CALL copy_fm_to_dbcsr(mo_notconv_fm,mo_notconv,error=error)!fm->dbcsr 00864 CALL cp_dbcsr_distribution_release (dist) 00865 CALL array_release (col_blk_size) 00866 CALL array_release (col_dist) 00867 00868 00869 CALL cp_dbcsr_multiply('n','n',1.0_dp,c_out,mo_notconv,0.0_dp,matrix_hc,& 00870 last_column=nmo_not_converged,error=error) 00871 CALL cp_dbcsr_multiply('n','n',1.0_dp,matrix_s,mo_notconv,0.0_dp,matrix_sc,& 00872 last_column=nmo_not_converged,error=error) 00873 00874 CALL cp_dbcsr_copy(matrix_z,matrix_sc,error=error) 00875 CALL cp_dbcsr_scale_by_vector(matrix_z,eig_not_conv,side='right',error=error) 00876 CALL cp_dbcsr_add(matrix_z,matrix_hc,-1.0_dp,1.0_dp,error=error) 00877 00878 DEALLOCATE(eig_not_conv,STAT=istat) 00879 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00880 00881 00882 ! matrix_mm 00883 CALL cp_create_bl_distribution (col_dist, col_blk_size, nmo_not_converged, & 00884 dbcsr_mp_npcols(dbcsr_distribution_mp(cp_dbcsr_distribution(matrix_s)))) 00885 CALL cp_create_bl_distribution (row_dist, row_blk_size, nmo_not_converged, & 00886 dbcsr_mp_nprows(dbcsr_distribution_mp(cp_dbcsr_distribution(matrix_s)))) 00887 CALL dbcsr_distribution_new (dist, dbcsr_distribution_mp(cp_dbcsr_distribution(matrix_s)),& 00888 row_dist,col_dist) 00889 CALL cp_dbcsr_init_p(matrix_mm,error=error) 00890 CALL cp_dbcsr_create(matrix_mm,"matrix_mm",dist,dbcsr_type_no_symmetry,& 00891 row_blk_size,col_blk_size,0,0,dbcsr_type_real_default,error=error) 00892 CALL cp_dbcsr_distribution_release (dist) 00893 CALL array_release (col_blk_size) 00894 CALL array_release (col_dist) 00895 CALL array_release (row_blk_size) 00896 CALL array_release (row_dist) 00897 00898 CALL cp_dbcsr_multiply('t','n',1.0_dp,mo_notconv,matrix_hc,0.0_dp,matrix_mm,& 00899 last_column=nmo_not_converged, error=error) 00900 00901 ELSE 00902 mo_notconv=>mo_coeff_b 00903 mo_notconv_fm=>mo_coeff 00904 c_out => matrix_h 00905 END IF 00906 00907 ! allocate matrix_pz using as template matrix_z 00908 CALL cp_dbcsr_init_p(matrix_pz,error=error) 00909 CALL cp_dbcsr_create(matrix_pz,"matrix_pz",cp_dbcsr_distribution(matrix_z),& 00910 dbcsr_type_no_symmetry,cp_dbcsr_row_block_sizes(matrix_z),& 00911 cp_dbcsr_col_block_sizes(matrix_z),0,0,dbcsr_type_real_default,error=error) 00912 00913 IF(do_apply_preconditioner) THEN 00914 IF(preconditioner%in_use/=0) THEN 00915 CALL apply_preconditioner(preconditioner,matrix_z,matrix_pz,error=error) 00916 ELSE 00917 CALL cp_dbcsr_copy(matrix_pz,matrix_z,error=error) 00918 END IF 00919 ELSE 00920 CALL cp_dbcsr_copy(matrix_pz,matrix_z,error=error) 00921 END IF 00922 00923 !allocate NMOxNMO full matrices 00924 nmat = nmo_not_converged 00925 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=nmat,ncol_global=nmat, & 00926 context=mo_coeff%matrix_struct%context, & 00927 para_env=mo_coeff%matrix_struct%para_env,error=error) 00928 CALL cp_fm_create(matrix_mm_fm,fm_struct_tmp,name="m_tmp_mxm",error=error) 00929 CALL cp_fm_create(matrix_mmt_fm,fm_struct_tmp,name="mt_tmp_mxm",error=error) 00930 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00931 00932 !allocate 2NMOx2NMO full matrices 00933 nmat2 = 2* nmo_not_converged 00934 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=nmat2,ncol_global=nmat2, & 00935 context=mo_coeff%matrix_struct%context, & 00936 para_env=mo_coeff%matrix_struct%para_env,error=error) 00937 00938 CALL cp_fm_create(s_block,fm_struct_tmp,name="sb",error=error) 00939 CALL cp_fm_create(h_block,fm_struct_tmp,name="hb",error=error) 00940 CALL cp_fm_create(v_block,fm_struct_tmp,name="vb",error=error) 00941 CALL cp_fm_create(w_block,fm_struct_tmp,name="wb",error=error) 00942 ALLOCATE(evals(nmat2)) 00943 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00944 00945 ! compute CSC 00946 CALL cp_fm_set_all(s_block,0.0_dp, 1.0_dp,error=error) 00947 ! compute CHC 00948 CALL copy_dbcsr_to_fm(matrix_mm,matrix_mm_fm,error=error) 00949 CALL cp_fm_to_fm_submat(matrix_mm_fm,h_block,nmat,nmat,1,1,1,1,error=error) 00950 00951 ! compute the bottom left ZSC (top right is transpose) 00952 CALL cp_dbcsr_multiply('t','n',1.0_dp,matrix_pz,matrix_sc,0.0_dp,matrix_mm,last_column=nmat,error=error) 00953 ! set the bottom left part of S[C,Z] block matrix ZSC 00954 !copy sparse to full 00955 CALL copy_dbcsr_to_fm(matrix_mm,matrix_mm_fm,error=error) 00956 CALL cp_fm_to_fm_submat(matrix_mm_fm,s_block,nmat,nmat,1,1,1+nmat,1,error=error) 00957 CALL cp_fm_transpose(matrix_mm_fm,matrix_mmt_fm,error=error) 00958 CALL cp_fm_to_fm_submat(matrix_mmt_fm,s_block,nmat,nmat,1,1,1,1+nmat,error=error) 00959 00960 ! compute the bottom left ZHC (top right is transpose) 00961 CALL cp_dbcsr_multiply('t','n',1.0_dp,matrix_pz,matrix_hc,0.0_dp,matrix_mm,last_column=nmat,error=error) 00962 ! set the bottom left part of S[C,Z] block matrix ZHC 00963 CALL copy_dbcsr_to_fm(matrix_mm,matrix_mm_fm,error=error) 00964 CALL cp_fm_to_fm_submat(matrix_mm_fm,h_block,nmat,nmat,1,1,1+nmat,1,error=error) 00965 CALL cp_fm_transpose(matrix_mm_fm,matrix_mmt_fm,error=error) 00966 CALL cp_fm_to_fm_submat(matrix_mmt_fm,h_block,nmat,nmat,1,1,1,1+nmat,error=error) 00967 00968 CALL cp_fm_release(matrix_mmt_fm,error=error) 00969 00970 ! (reuse matrix_sc and matrix_hc to computr HZ and SZ) 00971 CALL cp_dbcsr_get_info(matrix_pz,nfullrows_total=n,nfullcols_total=k) 00972 CALL cp_dbcsr_multiply('n','n',1.0_dp,c_out,matrix_pz,0.0_dp,matrix_hc,last_column=k,error=error) 00973 CALL cp_dbcsr_multiply('n','n',1.0_dp,matrix_s,matrix_pz,0.0_dp,matrix_sc,last_column=k,error=error) 00974 00975 ! compute the bottom right ZSZ 00976 CALL cp_dbcsr_multiply('t','n',1.0_dp,matrix_pz,matrix_sc,0.0_dp,matrix_mm,last_column=k,error=error) 00977 ! set the bottom right part of S[C,Z] block matrix ZSZ 00978 CALL copy_dbcsr_to_fm(matrix_mm,matrix_mm_fm,error=error) 00979 CALL cp_fm_to_fm_submat(matrix_mm_fm,s_block,nmat,nmat,1,1,1+nmat,1+nmat,error=error) 00980 00981 ! compute the bottom right ZHZ 00982 CALL cp_dbcsr_multiply('t','n',1.0_dp,matrix_pz,matrix_hc,0.0_dp,matrix_mm,last_column=k,error=error) 00983 ! set the bottom right part of H[C,Z] block matrix ZHZ 00984 CALL copy_dbcsr_to_fm(matrix_mm,matrix_mm_fm,error=error) 00985 CALL cp_fm_to_fm_submat(matrix_mm_fm,h_block,nmat,nmat,1,1,1+nmat,1+nmat,error=error) 00986 00987 CALL cp_dbcsr_release_p(matrix_mm, error=error) 00988 CALL cp_dbcsr_release_p(matrix_sc, error=error) 00989 CALL cp_dbcsr_release_p(matrix_hc, error=error) 00990 00991 CALL reduce_extended_space(s_block,h_block,v_block,w_block, evals, nmat2, error=error) 00992 00993 ! allocate two (nao x nmat) full matrix 00994 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=nao,ncol_global=nmat, & 00995 context=mo_coeff%matrix_struct%context, & 00996 para_env=mo_coeff%matrix_struct%para_env,error=error) 00997 CALL cp_fm_create(matrix_nm_fm,fm_struct_tmp,name="m_nxm",error=error) 00998 CALL cp_fm_create(matrix_z_fm,fm_struct_tmp,name="m_nxm",error=error) 00999 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 01000 01001 CALL copy_dbcsr_to_fm(matrix_pz,matrix_z_fm,error=error) 01002 ! extract egenvectors 01003 CALL cp_fm_to_fm_submat(v_block,matrix_mm_fm, nmat,nmat,1,1,1,1,error=error) 01004 CALL cp_fm_gemm('N','N',nao,nmat,nmat,1.0_dp,mo_notconv_fm,matrix_mm_fm,0.0_dp,matrix_nm_fm,error=error) 01005 CALL cp_fm_to_fm_submat(v_block,matrix_mm_fm, nmat,nmat,1+nmat,1,1,1,error=error) 01006 CALL cp_fm_gemm('N','N',nao,nmat,nmat,1.0_dp,matrix_z_fm,matrix_mm_fm,1.0_dp,matrix_nm_fm,error=error) 01007 01008 CALL cp_dbcsr_release_p(matrix_z, error=error) 01009 CALL cp_dbcsr_release_p(matrix_pz, error=error) 01010 CALL cp_fm_release(matrix_z_fm,error=error) 01011 CALL cp_fm_release(s_block,error=error) 01012 CALL cp_fm_release(h_block,error=error) 01013 CALL cp_fm_release(w_block,error=error) 01014 CALL cp_fm_release(v_block,error=error) 01015 CALL cp_fm_release(matrix_mm_fm,error=error) 01016 01017 01018 ! in case some vector are already converged only a subset of vectors are copied in the MOS 01019 IF (nmo_converged>0) THEN 01020 jj=1 01021 DO j=1,nset_not_conv 01022 i_first=inotconv_set(j,1) 01023 i_last=inotconv_set(j,2) 01024 n=i_last-i_first+1 01025 CALL cp_fm_to_fm_submat(matrix_nm_fm,mo_coeff,nao,n,1,jj,1,i_first,error=error) 01026 eigenvalues(i_first:i_last) = evals(jj:jj+n-1) 01027 jj=jj+n 01028 END DO 01029 DEALLOCATE(iconv_set,STAT=istat) 01030 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01031 DEALLOCATE (inotconv_set,STAT=istat) 01032 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01033 01034 CALL cp_dbcsr_release_p(mo_notconv, error=error) 01035 CALL cp_dbcsr_release_p(c_out, error=error) 01036 CALL cp_fm_release(mo_notconv_fm,error=error) 01037 ELSE 01038 CALL cp_fm_to_fm(matrix_nm_fm,mo_coeff,error=error) 01039 eigenvalues(1:nmo) = evals(1:nmo) 01040 END IF 01041 DEALLOCATE(evals,STAT=istat) 01042 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01043 01044 CALL cp_fm_release(matrix_nm_fm,error=error) 01045 CALL copy_fm_to_dbcsr(mo_coeff,mo_coeff_b,error=error)!fm->dbcsr 01046 01047 t2 = m_walltime() 01048 IF (output_unit > 0) THEN 01049 WRITE(output_unit,'(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)')& 01050 iter, nmo_converged, max_norm, min_norm, t2-t1 01051 END IF 01052 t1=m_walltime() 01053 01054 01055 END DO ! iter 01056 01057 01058 DEALLOCATE(ritz_coeff,STAT=istat) 01059 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01060 DEALLOCATE(iconv,STAT=istat) 01061 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01062 DEALLOCATE(inotconv,STAT=istat) 01063 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01064 DEALLOCATE(vnorm,STAT=istat) 01065 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01066 01067 CALL timestop(handle) 01068 01069 END SUBROUTINE generate_extended_space_sparse 01070 01071 ! ***************************************************************************** 01072 ! ***************************************************************************** 01073 01074 SUBROUTINE reduce_extended_space(s_block,h_block,v_block,w_block, evals, ndim, error) 01075 01076 TYPE(cp_fm_type), POINTER :: s_block, h_block, v_block, 01077 w_block 01078 REAL(dp), DIMENSION(:) :: evals 01079 INTEGER :: ndim 01080 TYPE(cp_error_type), INTENT(inout) :: error 01081 01082 CHARACTER(len=*), PARAMETER :: routineN = 'reduce_extended_space', 01083 routineP = moduleN//':'//routineN 01084 01085 INTEGER :: handle, info 01086 LOGICAL :: failure 01087 01088 failure=.FALSE. 01089 01090 CALL timeset(routineN,handle) 01091 01092 CALL cp_fm_to_fm(s_block,w_block,error=error) 01093 CALL cp_fm_cholesky_decompose(s_block,info_out=info,error=error) 01094 IF(info==0) THEN 01095 CALL cp_fm_triangular_invert(s_block,error=error) 01096 CALL cp_fm_cholesky_restore(H_block,ndim,S_block,w_block,"MULTIPLY",pos="RIGHT",error=error) 01097 CALL cp_fm_cholesky_restore(w_block,ndim,S_block,H_block,"MULTIPLY",pos="LEFT",transa="T",error=error) 01098 CALL choose_eigv_solver(H_block,w_block,evals,error=error) 01099 CALL cp_fm_cholesky_restore(w_block,ndim,S_block,v_block,"MULTIPLY",error=error) 01100 ELSE 01101 ! S^(-1/2) 01102 CALL cp_fm_power(w_block,s_block,-0.5_dp,1.0E-5_dp,info,error=error) 01103 CALL cp_fm_to_fm(w_block,s_block,error=error) 01104 CALL cp_fm_gemm('N','N',ndim,ndim,ndim,1.0_dp,H_block,s_block,0.0_dp,w_block,error=error) 01105 CALL cp_fm_gemm('N','N',ndim,ndim,ndim,1.0_dp,s_block,w_block,0.0_dp,H_block,error=error) 01106 CALL choose_eigv_solver(H_block,w_block,evals,error=error) 01107 CALL cp_fm_gemm('N','N',ndim,ndim,ndim,1.0_dp,s_block,w_block,0.0_dp,v_block,error=error) 01108 END IF 01109 01110 CALL timestop(handle) 01111 01112 END SUBROUTINE reduce_extended_space 01113 01114 01115 END MODULE qs_scf_block_davidson
1.7.3