|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00013 MODULE qs_ot 00014 USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,& 00015 cp_dbcsr_cholesky_invert,& 00016 cp_dbcsr_cholesky_restore 00017 USE cp_dbcsr_diag, ONLY: cp_dbcsr_heevd,& 00018 cp_dbcsr_syevd 00019 USE cp_dbcsr_interface, ONLY: & 00020 cp_dbcsr_add, cp_dbcsr_add_on_diag, cp_dbcsr_copy, & 00021 cp_dbcsr_distribution, cp_dbcsr_filter, cp_dbcsr_frobenius_norm, & 00022 cp_dbcsr_gershgorin_norm, cp_dbcsr_get_block_p, cp_dbcsr_get_info, & 00023 cp_dbcsr_get_occupation, cp_dbcsr_hadamard_product, cp_dbcsr_init, & 00024 cp_dbcsr_init_p, cp_dbcsr_iterator_blocks_left, & 00025 cp_dbcsr_iterator_next_block, cp_dbcsr_iterator_start, & 00026 cp_dbcsr_iterator_stop, cp_dbcsr_multiply, cp_dbcsr_release, & 00027 cp_dbcsr_release_p, cp_dbcsr_scale, cp_dbcsr_scale_by_vector, & 00028 cp_dbcsr_set, cp_dbcsr_transposed 00029 USE cp_dbcsr_types, ONLY: cp_dbcsr_iterator,& 00030 cp_dbcsr_type 00031 USE dbcsr_methods, ONLY: dbcsr_distribution_mp,& 00032 dbcsr_mp_group 00033 USE f77_blas 00034 USE kinds, ONLY: dp 00035 USE message_passing, ONLY: mp_sum 00036 USE preconditioner, ONLY: apply_preconditioner 00037 USE preconditioner_types, ONLY: preconditioner_type 00038 USE qs_ot_types, ONLY: qs_ot_type 00039 USE scp_coeff_types, ONLY: aux_coeff_set_type,& 00040 aux_coeff_type,& 00041 get_aux_coeff 00042 USE termination, ONLY: stop_program 00043 USE timings, ONLY: timeset,& 00044 timestop 00045 #include "cp_common_uses.h" 00046 00047 IMPLICIT NONE 00048 PRIVATE 00049 00050 PUBLIC :: qs_ot_get_p 00051 PUBLIC :: qs_ot_get_orbitals 00052 PUBLIC :: qs_ot_get_derivative 00053 PUBLIC :: qs_ot_get_orbitals_ref 00054 PUBLIC :: qs_ot_get_derivative_ref 00055 PUBLIC :: qs_ot_new_preconditioner 00056 PUBLIC :: qs_ot_get_scp_dft_derivative 00057 PUBLIC :: qs_ot_get_scp_dft_coeffs 00058 PUBLIC :: qs_ot_get_scp_nddo_derivative 00059 PUBLIC :: qs_ot_get_scp_nddo_coeffs 00060 PRIVATE :: qs_ot_p2m_diag 00061 PRIVATE :: qs_ot_sinc 00062 PRIVATE :: qs_ot_ref_poly 00063 PRIVATE :: qs_ot_ref_chol 00064 PRIVATE :: qs_ot_ref_lwdn 00065 PRIVATE :: qs_ot_ref_decide 00066 PRIVATE :: qs_ot_ref_update 00067 PRIVATE :: qs_ot_refine 00068 PRIVATE :: qs_ot_on_the_fly_localize 00069 00070 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ot' 00071 00072 CONTAINS 00073 00074 ! gets ready to use the preconditioner/ or renew the preconditioner 00075 ! only keeps a pointer to the preconditioner. 00076 ! If you change the preconditioner, you have to call this routine 00077 ! you remain responsible of proper deallocate of your preconditioner 00078 ! (or you can reuse it on the next step of the computation) 00079 ! ***************************************************************************** 00080 SUBROUTINE qs_ot_new_preconditioner(qs_ot_env,preconditioner,error) 00081 TYPE(qs_ot_type) :: qs_ot_env 00082 TYPE(preconditioner_type), POINTER :: preconditioner 00083 TYPE(cp_error_type), INTENT(inout) :: error 00084 00085 CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_new_preconditioner', 00086 routineP = moduleN//':'//routineN 00087 00088 INTEGER :: istat, ncoef 00089 LOGICAL :: failure, mixed_precision 00090 00091 failure = .FALSE. 00092 00093 qs_ot_env%preconditioner => preconditioner 00094 qs_ot_env%os_valid = .FALSE. 00095 IF (.NOT. ASSOCIATED(qs_ot_env%matrix_psc0)) THEN 00096 CALL cp_dbcsr_init_p(qs_ot_env%matrix_psc0, error=error) 00097 CALL cp_dbcsr_copy(qs_ot_env%matrix_psc0,qs_ot_env%matrix_sc0,'matrix_psc0',error=error) 00098 ENDIF 00099 00100 mixed_precision = qs_ot_env%settings%mixed_precision 00101 00102 IF (.NOT. qs_ot_env%use_dx) THEN 00103 qs_ot_env%use_dx=.TRUE. 00104 CALL cp_dbcsr_init_p(qs_ot_env%matrix_dx, error=error) 00105 CALL cp_dbcsr_copy(qs_ot_env%matrix_dx,qs_ot_env%matrix_gx,'matrix_dx',error=error) 00106 IF (qs_ot_env%settings%do_rotation) THEN 00107 CALL cp_dbcsr_init_p(qs_ot_env%rot_mat_dx, error=error) 00108 CALL cp_dbcsr_copy(qs_ot_env%rot_mat_dx,qs_ot_env%rot_mat_gx,'rot_mat_dx',error=error) 00109 ENDIF 00110 IF (qs_ot_env%settings%do_ener) THEN 00111 ncoef = SIZE ( qs_ot_env % ener_gx) 00112 ALLOCATE ( qs_ot_env%ener_dx ( ncoef ), STAT = istat ) 00113 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00114 qs_ot_env%ener_dx = 0.0_dp 00115 ENDIF 00116 ! ***SCP 00117 IF ( qs_ot_env % settings % scp_dft ) THEN 00118 ncoef = SIZE ( qs_ot_env % gx ) 00119 ALLOCATE ( qs_ot_env%dx ( ncoef ), STAT = istat ) 00120 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00121 qs_ot_env%dx = 0.0_dp 00122 ENDIF 00123 IF ( qs_ot_env % settings % scp_nddo ) THEN 00124 ALLOCATE(qs_ot_env % dxmat) 00125 CALL cp_dbcsr_init(qs_ot_env % dxmat, error=error) 00126 CALL cp_dbcsr_copy(qs_ot_env % dxmat, qs_ot_env % gxmat, "SCP_DXMAT", error=error) 00127 CALL cp_dbcsr_set (qs_ot_env % dxmat, 0.0_dp, error=error ) 00128 ENDIF 00129 ! ***SCP 00130 ENDIF 00131 00132 END SUBROUTINE qs_ot_new_preconditioner 00133 00134 ! ***************************************************************************** 00135 SUBROUTINE qs_ot_on_the_fly_localize(qs_ot_env, C_NEW, SC, G_OLD, D, error) 00136 ! 00137 TYPE(qs_ot_type) :: qs_ot_env 00138 TYPE(cp_dbcsr_type), POINTER :: C_NEW, SC, G_OLD, D 00139 TYPE(cp_error_type), INTENT(INOUT) :: error 00140 00141 CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_on_the_fly_localize', 00142 routineP = moduleN//':'//routineN 00143 INTEGER, PARAMETER :: taylor_order = 50 00144 REAL(KIND=dp), PARAMETER :: alpha = 0.1_dp, 00145 f2_eps = 0.01_dp, 00146 rone = 1.0_dp, rzero = 0.0_dp 00147 00148 INTEGER :: blk, col, col_size, handle, 00149 i, k, n, output_unit, p, row, 00150 row_size 00151 REAL(dp), DIMENSION(:, :), POINTER :: DATA 00152 REAL(KIND=dp) :: expfactor, f2, norm_fro, 00153 norm_gct, tmp 00154 TYPE(cp_dbcsr_iterator) :: iter 00155 TYPE(cp_dbcsr_type), POINTER :: C, Gp1, Gp2, GU, U 00156 TYPE(cp_logger_type), POINTER :: logger 00157 00158 CALL timeset(routineN,handle) 00159 ! 00160 ! 00161 CALL cp_dbcsr_get_info(C_NEW,nfullrows_total=n,nfullcols_total=k) 00162 ! 00163 ! C = C*expm(-G) 00164 GU => qs_ot_env%buf1_k_k_nosym ! a buffer 00165 U => qs_ot_env%buf2_k_k_nosym ! a buffer 00166 Gp1=> qs_ot_env%buf3_k_k_nosym ! a buffer 00167 Gp2=> qs_ot_env%buf4_k_k_nosym ! a buffer 00168 C => qs_ot_env%buf1_n_k ! a buffer 00169 ! 00170 ! compute the derivative of the norm 00171 !------------------------------------------------------------------- 00172 ! (x^2+eps)^1/2 00173 f2 = 0.0_dp 00174 CALL cp_dbcsr_copy(C,C_NEW,error=error) 00175 CALL cp_dbcsr_iterator_start(iter, C) 00176 DO WHILE (cp_dbcsr_iterator_blocks_left (iter)) 00177 CALL cp_dbcsr_iterator_next_block(iter, row, col, DATA, blk,& 00178 row_size=row_size, col_size=col_size) 00179 DO p=1,col_size! p 00180 DO i=1,row_size! i 00181 tmp = SQRT( DATA(i,p)**2 + f2_eps ) 00182 f2 = f2 + tmp 00183 DATA(i,p) = DATA(i,p) / tmp 00184 ENDDO 00185 ENDDO 00186 ENDDO 00187 CALL cp_dbcsr_iterator_stop(iter) 00188 CALL mp_sum(f2,dbcsr_mp_group(dbcsr_distribution_mp(cp_dbcsr_distribution(C)))) 00189 logger => cp_error_get_logger(error) 00190 output_unit = cp_logger_get_default_io_unit(logger) 00191 IF(output_unit>0) WRITE(output_unit,*) routineN//' f2 =',f2 00192 ! 00193 ! 00194 CALL cp_dbcsr_multiply('T','N',1.0_dp,C,C_NEW,0.0_dp,GU,error=error) 00195 ! 00196 ! antisymetrize 00197 CALL cp_dbcsr_transposed (U, GU, shallow_data_copy=.FALSE., & 00198 use_distribution=cp_dbcsr_distribution(GU), & 00199 transpose_distribution=.FALSE., & 00200 error=error) 00201 CALL cp_dbcsr_add(GU,U,alpha_scalar=-0.5_dp,beta_scalar=0.5_dp,error=error) 00202 !------------------------------------------------------------------- 00203 ! 00204 norm_fro = cp_dbcsr_frobenius_norm(GU) 00205 norm_gct = cp_dbcsr_gershgorin_norm(GU) 00206 !write(*,*) 'qs_ot_localize: ||P-I||_f=',norm_fro,' ||P-I||_GCT=',norm_gct 00207 ! 00208 !kscale = CEILING(LOG(MIN(norm_fro,norm_gct))/LOG(2.0_dp)) 00209 !scale = LOG(MIN(norm_fro,norm_gct))/LOG(2.0_dp) 00210 !write(*,*) 'qs_ot_localize: scale=',scale,' kscale=',kscale 00211 ! 00212 ! rescale for steepest descent 00213 CALL cp_dbcsr_scale(GU, -alpha, error=error) 00214 ! 00215 ! compute unitary transform 00216 ! zeroth and first order 00217 expfactor = 1.0_dp 00218 CALL cp_dbcsr_copy(U,GU,error=error) 00219 CALL cp_dbcsr_scale(U,expfactor,error=error) 00220 CALL cp_dbcsr_add_on_diag(U,1.0_dp,error=error) 00221 ! other orders 00222 CALL cp_dbcsr_copy(Gp1,GU,error=error) 00223 DO i = 2,taylor_order 00224 ! new power of G 00225 CALL cp_dbcsr_multiply('N','N',1.0_dp,GU,Gp1,0.0_dp,Gp2,error=error) 00226 CALL cp_dbcsr_copy(Gp1,Gp2,error=error) 00227 ! add to the taylor expansion so far 00228 expfactor = expfactor / REAL(i,KIND=dp) 00229 CALL cp_dbcsr_add(U,Gp1,alpha_scalar=1.0_dp,beta_scalar=expfactor,error=error) 00230 norm_fro = cp_dbcsr_frobenius_norm(Gp1) 00231 !write(*,*) 'Taylor expansion i=',i,' norm(X^i)/i!=',norm_fro*expfactor 00232 IF(norm_fro*expfactor.LT.1.0E-10_dp) EXIT 00233 ENDDO 00234 ! 00235 ! rotate MOs 00236 CALL cp_dbcsr_multiply('N','N',1.0_dp,C_NEW,U,0.0_dp,C,error=error) 00237 CALL cp_dbcsr_copy(C_NEW,C,error=error) 00238 ! 00239 ! rotate SC 00240 CALL cp_dbcsr_multiply('N','N',1.0_dp,SC,U,0.0_dp,C,error=error) 00241 CALL cp_dbcsr_copy(SC,C,error=error) 00242 ! 00243 ! rotate D_i 00244 CALL cp_dbcsr_multiply('N','N',1.0_dp,D,U,0.0_dp,C,error=error) 00245 CALL cp_dbcsr_copy(D,C,error=error) 00246 ! 00247 ! rotate G_i-1 00248 IF(ASSOCIATED(G_OLD)) THEN 00249 CALL cp_dbcsr_multiply('N','N',1.0_dp,G_OLD,U,0.0_dp,C,error=error) 00250 CALL cp_dbcsr_copy(G_OLD,C,error=error) 00251 ENDIF 00252 ! 00253 CALL timestop(handle) 00254 END SUBROUTINE qs_ot_on_the_fly_localize 00255 00256 ! ***************************************************************************** 00257 SUBROUTINE qs_ot_ref_chol(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, update, error) 00258 ! 00259 TYPE(qs_ot_type) :: qs_ot_env 00260 TYPE(cp_dbcsr_type) :: C_OLD, C_TMP, C_NEW, P, SC 00261 LOGICAL, INTENT(IN) :: update 00262 TYPE(cp_error_type), INTENT(INOUT) :: error 00263 00264 CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_ref_chol', 00265 routineP = moduleN//':'//routineN 00266 00267 INTEGER :: handle, k, n 00268 00269 CALL timeset(routineN,handle) 00270 ! 00271 CALL cp_dbcsr_get_info(C_NEW,nfullrows_total=n,nfullcols_total=k) 00272 ! 00273 ! P = U'*U 00274 CALL cp_dbcsr_cholesky_decompose(P,k,qs_ot_env%para_env,qs_ot_env%blacs_env,error=error) 00275 ! 00276 ! C_NEW = C_OLD*inv(U) 00277 CALL cp_dbcsr_cholesky_restore(C_OLD,k,P,C_NEW,op="SOLVE",pos="RIGHT",& 00278 transa="N",para_env=qs_ot_env%para_env,blacs_env=qs_ot_env%blacs_env,& 00279 error=error) 00280 ! 00281 ! Update SC if needed 00282 IF(update) THEN 00283 CALL cp_dbcsr_cholesky_restore(SC,k,P,C_TMP,op="SOLVE",pos="RIGHT",& 00284 transa="N",para_env=qs_ot_env%para_env,blacs_env=qs_ot_env%blacs_env,error=error) 00285 CALL cp_dbcsr_copy(SC,C_TMP,error=error) 00286 ENDIF 00287 ! 00288 CALL timestop(handle) 00289 END SUBROUTINE qs_ot_ref_chol 00290 00291 ! ***************************************************************************** 00292 SUBROUTINE qs_ot_ref_lwdn(qs_ot_env, C_OLD, C_TMP, C_NEW, P, SC, update, error) 00293 ! 00294 TYPE(qs_ot_type) :: qs_ot_env 00295 TYPE(cp_dbcsr_type) :: C_OLD, C_TMP, C_NEW, P, SC 00296 LOGICAL, INTENT(IN) :: update 00297 TYPE(cp_error_type), INTENT(INOUT) :: error 00298 00299 CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_ref_lwdn', 00300 routineP = moduleN//':'//routineN 00301 00302 INTEGER :: handle, i, istat, k, n 00303 LOGICAL :: failure 00304 REAL(dp), ALLOCATABLE, DIMENSION(:) :: eig, fun 00305 TYPE(cp_dbcsr_type), POINTER :: V, W 00306 00307 failure = .FALSE. 00308 00309 CALL timeset(routineN,handle) 00310 ! 00311 CALL cp_dbcsr_get_info(C_NEW,nfullrows_total=n,nfullcols_total=k) 00312 ! 00313 V => qs_ot_env%buf1_k_k_nosym ! a buffer 00314 W => qs_ot_env%buf2_k_k_nosym ! a buffer 00315 ALLOCATE(eig(k), fun(k), STAT=istat) 00316 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00317 ! 00318 CALL cp_dbcsr_syevd(P,V,eig,qs_ot_env%para_env,qs_ot_env%blacs_env,error=error) 00319 ! 00320 ! compute the P^(-1/2) 00321 DO i = 1,k 00322 IF(eig(i).LE.0.0_dp) & 00323 & CALL stop_program(routineN,moduleN,__LINE__,"P not positive definite") 00324 IF(eig(i).LT.1.0E-8_dp) THEN 00325 fun(i)=0.0_dp 00326 ELSE 00327 fun(i)=1.0_dp/SQRT(eig(i)) 00328 ENDIF 00329 ENDDO 00330 CALL cp_dbcsr_copy(W,V,error=error) 00331 CALL cp_dbcsr_scale_by_vector(V,alpha=fun,side='right',error=error) 00332 CALL cp_dbcsr_multiply('N','T',1.0_dp,W,V,0.0_dp,P,error=error) 00333 ! 00334 ! Update C 00335 CALL cp_dbcsr_multiply('N','N',1.0_dp,C_OLD,P,0.0_dp,C_NEW,error=error) 00336 ! 00337 ! Update SC if needed 00338 IF(update) THEN 00339 CALL cp_dbcsr_multiply('N','N',1.0_dp,SC,P,0.0_dp,C_TMP,error=error) 00340 CALL cp_dbcsr_copy(SC,C_TMP,error=error) 00341 ENDIF 00342 ! 00343 DEALLOCATE(eig, fun, STAT=istat) 00344 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00345 ! 00346 CALL timestop(handle) 00347 END SUBROUTINE qs_ot_ref_lwdn 00348 00349 ! ***************************************************************************** 00350 SUBROUTINE qs_ot_ref_poly(qs_ot_env,C_OLD,C_TMP,C_NEW,P,SC,norm_in,update,output_unit,error) 00351 ! 00352 TYPE(qs_ot_type) :: qs_ot_env 00353 TYPE(cp_dbcsr_type), POINTER :: C_OLD, C_TMP, C_NEW, P 00354 TYPE(cp_dbcsr_type) :: SC 00355 REAL(dp), INTENT(IN) :: norm_in 00356 LOGICAL, INTENT(IN) :: update 00357 INTEGER, INTENT(IN) :: output_unit 00358 TYPE(cp_error_type), INTENT(INOUT) :: error 00359 00360 CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_ref_poly', 00361 routineP = moduleN//':'//routineN 00362 00363 INTEGER :: handle, irefine, k, n 00364 LOGICAL :: quick_exit 00365 REAL(dp) :: norm, norm_fro, norm_gct, 00366 occ_in, occ_out, rescale 00367 TYPE(cp_dbcsr_type), POINTER :: BUF1, BUF2, BUF_NOSYM, FT, FY 00368 00369 CALL timeset(routineN,handle) 00370 ! 00371 CALL cp_dbcsr_get_info(C_NEW,nfullrows_total=n,nfullcols_total=k) 00372 ! 00373 BUF_NOSYM => qs_ot_env%buf1_k_k_nosym! a buffer 00374 BUF1 => qs_ot_env%buf1_k_k_sym ! a buffer 00375 BUF2 => qs_ot_env%buf2_k_k_sym ! a buffer 00376 FY => qs_ot_env%buf3_k_k_sym ! a buffer 00377 FT => qs_ot_env%buf4_k_k_sym ! a buffer 00378 ! 00379 ! initialize the norm (already computed in qs_ot_get_orbitals_ref) 00380 norm = norm_in 00381 ! 00382 ! can we do a quick exit? 00383 quick_exit = .FALSE. 00384 IF(norm.LT.qs_ot_env%settings%eps_irac_quick_exit) quick_exit = .TRUE. 00385 ! 00386 ! lets refine 00387 rescale = 1.0_dp 00388 DO irefine = 1,qs_ot_env%settings%max_irac 00389 ! 00390 ! rescaling 00391 IF(norm.GT.1.0_dp) THEN 00392 IF(output_unit>0) WRITE(output_unit,'(A,I3,A)') & 00393 & routineN,irefine,': we rescale (C+a*D)' 00394 CALL cp_dbcsr_scale(P,1.0_dp/norm,error=error) 00395 rescale = rescale/SQRT(norm) 00396 ENDIF 00397 ! 00398 ! get the refinement polynomial 00399 CALL qs_ot_refine(P, FY, BUF1, BUF2, qs_ot_env%settings%irac_degree, & 00400 qs_ot_env%settings%eps_irac_filter_matrix, output_unit, error) 00401 ! 00402 ! collect the transformation 00403 IF(irefine.EQ.1) THEN 00404 CALL cp_dbcsr_copy(FT, FY, name='FT', error=error) 00405 ELSE 00406 CALL cp_dbcsr_multiply('N', 'N', 1.0_dp, FT, FY, 0.0_dp, BUF1, error=error) 00407 IF (qs_ot_env%settings%eps_irac_filter_matrix.GT.0.0_dp) THEN 00408 occ_in = cp_dbcsr_get_occupation(buf1) 00409 CALL cp_dbcsr_filter(buf1,qs_ot_env%settings%eps_irac_filter_matrix,error=error) 00410 occ_out = cp_dbcsr_get_occupation(buf1) 00411 IF(output_unit>0.AND.qs_ot_env%settings%eps_irac_filter_matrix.GT.0.0_dp) & 00412 WRITE(output_unit,'(2(A,F8.5))') routinen//' filter(BUF1): occ_in',occ_in,' occ_out',occ_out 00413 ENDIF 00414 CALL cp_dbcsr_copy(FT, BUF1, name='FT', error=error) 00415 ENDIF 00416 ! 00417 ! quick exit if possible 00418 IF(quick_exit) THEN 00419 IF(output_unit>0) WRITE(output_unit,'(A,I3,A)') & 00420 & routineN,irefine,': quick exit!' 00421 EXIT 00422 ENDIF 00423 ! 00424 ! P = FY^T * P * FY 00425 CALL cp_dbcsr_multiply('N', 'N', 1.0_dp, P, FY, 0.0_dp, BUF_NOSYM, error=error) 00426 IF (qs_ot_env%settings%eps_irac_filter_matrix.GT.0.0_dp) THEN 00427 occ_in = cp_dbcsr_get_occupation(buf_nosym) 00428 CALL cp_dbcsr_filter(buf_nosym,qs_ot_env%settings%eps_irac_filter_matrix,error=error) 00429 occ_out = cp_dbcsr_get_occupation(buf_nosym) 00430 IF(output_unit>0.AND.qs_ot_env%settings%eps_irac_filter_matrix.GT.0.0_dp) & 00431 WRITE(output_unit,'(2(A,F8.5))') routinen//' filter(BUF_NOSYM): occ_in',occ_in,' occ_out',occ_out 00432 ENDIF 00433 CALL cp_dbcsr_multiply('N', 'N', 1.0_dp, FY, BUF_NOSYM, 0.0_dp, P, error=error) 00434 IF (qs_ot_env%settings%eps_irac_filter_matrix.GT.0.0_dp) THEN 00435 occ_in = cp_dbcsr_get_occupation(p) 00436 CALL cp_dbcsr_filter(p,qs_ot_env%settings%eps_irac_filter_matrix,error=error) 00437 occ_out = cp_dbcsr_get_occupation(p) 00438 IF(output_unit>0.AND.qs_ot_env%settings%eps_irac_filter_matrix.GT.0.0_dp) & 00439 WRITE(output_unit,'(2(A,F8.5))') routinen//' filter(P): occ_in',occ_in,' occ_out',occ_out 00440 ENDIF 00441 ! 00442 ! check ||P-1||_gct 00443 CALL cp_dbcsr_add_on_diag(P, -1.0_dp, error=error) 00444 norm_fro = cp_dbcsr_frobenius_norm(P) 00445 norm_gct = cp_dbcsr_gershgorin_norm(P) 00446 CALL cp_dbcsr_add_on_diag(P, 1.0_dp, error=error) 00447 norm = MIN(norm_gct,norm_fro) 00448 ! 00449 ! printing 00450 IF(output_unit>0) WRITE(output_unit,'(A,I3,A,E12.5)') & 00451 & routineN,irefine,': ||P-I||=',norm 00452 ! 00453 ! blows up 00454 IF (norm > 1.0E10_dp) THEN 00455 CALL stop_program(routineN,moduleN,__LINE__,& 00456 "Refinement blows up! "//& 00457 "We need you to improve the code, please post your input on "//& 00458 "the forum http://www.cp2k.org/") 00459 END IF 00460 ! 00461 ! can we do a quick exit next step? 00462 IF(norm.LT.qs_ot_env%settings%eps_irac_quick_exit) quick_exit = .TRUE. 00463 ! 00464 ! are we done? 00465 IF(norm.LT.qs_ot_env%settings%eps_irac) EXIT 00466 ! 00467 ENDDO 00468 ! 00469 ! C_NEW = C_NEW * FT * rescale 00470 CALL cp_dbcsr_multiply('N', 'N', rescale, C_OLD, FT, 0.0_dp, C_NEW, error=error) 00471 IF (qs_ot_env%settings%eps_irac_filter_matrix.GT.0.0_dp) THEN 00472 occ_in = cp_dbcsr_get_occupation(c_new) 00473 CALL cp_dbcsr_filter(c_new,qs_ot_env%settings%eps_irac_filter_matrix,error=error) 00474 occ_out = cp_dbcsr_get_occupation(c_new) 00475 IF(output_unit>0.AND.qs_ot_env%settings%eps_irac_filter_matrix.GT.0.0_dp) & 00476 WRITE(output_unit,'(2(A,F8.5))') routinen//' filter(C_NEW): occ_in',occ_in,' occ_out',occ_out 00477 ENDIF 00478 ! 00479 ! update SC = SC * FY * rescale 00480 IF(update) THEN 00481 CALL cp_dbcsr_multiply('N', 'N', rescale, SC, FT, 0.0_dp, C_TMP, error=error) 00482 IF (qs_ot_env%settings%eps_irac_filter_matrix.GT.0.0_dp) THEN 00483 occ_in = cp_dbcsr_get_occupation(c_tmp) 00484 CALL cp_dbcsr_filter(c_tmp,qs_ot_env%settings%eps_irac_filter_matrix,error=error) 00485 occ_out = cp_dbcsr_get_occupation(c_tmp) 00486 IF(output_unit>0.AND.qs_ot_env%settings%eps_irac_filter_matrix.GT.0.0_dp) & 00487 WRITE(output_unit,'(2(A,F8.5))') routinen//' filter(C_TMP): occ_in',occ_in,' occ_out',occ_out 00488 ENDIF 00489 CALL cp_dbcsr_copy(SC, C_TMP, error=error) 00490 ENDIF 00491 ! 00492 CALL timestop(handle) 00493 END SUBROUTINE qs_ot_ref_poly 00494 00495 ! ***************************************************************************** 00496 FUNCTION qs_ot_ref_update(qs_ot_env1) RESULT(update) 00497 ! 00498 TYPE(qs_ot_type) :: qs_ot_env1 00499 LOGICAL :: update 00500 00501 CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_ref_update', 00502 routineP = moduleN//':'//routineN 00503 00504 update = .FALSE. 00505 SELECT CASE(qs_ot_env1%settings%ot_method) 00506 CASE("CG") 00507 SELECT CASE(qs_ot_env1%settings%line_search_method) 00508 CASE("2PNT") 00509 IF(qs_ot_env1%line_search_count.EQ.2) update = .TRUE. 00510 CASE DEFAULT 00511 CALL stop_program(routineN,moduleN,__LINE__,"NYI") 00512 END SELECT 00513 CASE("DIIS") 00514 update = .TRUE. 00515 CASE DEFAULT 00516 CALL stop_program(routineN,moduleN,__LINE__,"NYI") 00517 END SELECT 00518 END FUNCTION qs_ot_ref_update 00519 00520 ! ***************************************************************************** 00521 SUBROUTINE qs_ot_ref_decide(qs_ot_env1, norm_in, ortho_irac) 00522 ! 00523 TYPE(qs_ot_type) :: qs_ot_env1 00524 REAL(dp), INTENT(IN) :: norm_in 00525 CHARACTER(LEN=*), INTENT(INOUT) :: ortho_irac 00526 00527 ortho_irac = qs_ot_env1%settings%ortho_irac 00528 IF(norm_in.LT.qs_ot_env1%settings%eps_irac_switch) ortho_irac = "POLY" 00529 END SUBROUTINE qs_ot_ref_decide 00530 00531 ! ***************************************************************************** 00532 SUBROUTINE qs_ot_get_orbitals_ref(matrix_c, matrix_s, matrix_x, matrix_sx, & 00533 & matrix_gx_old, matrix_dx, qs_ot_env, & 00534 & qs_ot_env1, output_unit, error) 00535 ! 00536 TYPE(cp_dbcsr_type), POINTER :: matrix_c, matrix_s, matrix_x, 00537 matrix_sx, matrix_gx_old, 00538 matrix_dx 00539 TYPE(qs_ot_type) :: qs_ot_env, qs_ot_env1 00540 INTEGER, INTENT(IN) :: output_unit 00541 TYPE(cp_error_type), INTENT(INOUT) :: error 00542 00543 CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_orbitals_ref', 00544 routineP = moduleN//':'//routineN 00545 REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp 00546 00547 CHARACTER(LEN=4) :: ortho_irac 00548 INTEGER :: handle, k, n 00549 LOGICAL :: on_the_fly_loc, update 00550 REAL(dp) :: norm, norm_fro, norm_gct, 00551 occ_in, occ_out 00552 TYPE(cp_dbcsr_type), POINTER :: C_NEW, C_OLD, C_TMP, D, 00553 G_OLD, P, S, SC 00554 00555 CALL timeset(routineN,handle) 00556 00557 IF(output_unit>0.AND.qs_ot_env1%settings%eps_irac_filter_matrix.GT.0.0_dp) & 00558 WRITE(output_unit,*) routinen//' eps_irac_filter_matrix=',& 00559 qs_ot_env1%settings%eps_irac_filter_matrix 00560 00561 CALL cp_dbcsr_get_info(matrix_c,nfullrows_total=n,nfullcols_total=k) 00562 ! 00563 C_NEW => matrix_c 00564 C_OLD => matrix_x ! need to be carefully updated for the gradient ! 00565 SC => matrix_sx! need to be carefully updated for the gradient ! 00566 G_OLD => matrix_gx_old ! need to be carefully updated for localization ! 00567 D => matrix_dx ! need to be carefully updated for localization ! 00568 S => matrix_s 00569 00570 P => qs_ot_env%p_k_k_sym ! a buffer 00571 C_TMP => qs_ot_env%buf1_n_k ! a buffer 00572 ! 00573 ! do we need to update C_OLD and SC? 00574 update = qs_ot_ref_update(qs_ot_env1) 00575 ! 00576 ! do we want to on the fly localize? 00577 ! for the moment this is set from the input, 00578 ! later we might want to localize every n-step or 00579 ! when the sparsity increases... 00580 on_the_fly_loc = qs_ot_env1%settings%on_the_fly_loc 00581 ! 00582 ! compute SC = S*C 00583 IF(ASSOCIATED(S)) THEN 00584 CALL cp_dbcsr_multiply('N','N',1.0_dp,S,C_OLD,0.0_dp,SC,error=error) 00585 IF (qs_ot_env1%settings%eps_irac_filter_matrix.GT.0.0_dp) THEN 00586 occ_in = cp_dbcsr_get_occupation(sc) 00587 CALL cp_dbcsr_filter(sc,qs_ot_env1%settings%eps_irac_filter_matrix,error=error) 00588 occ_out = cp_dbcsr_get_occupation(sc) 00589 IF(output_unit>0.AND.qs_ot_env1%settings%eps_irac_filter_matrix.GT.0.0_dp) & 00590 WRITE(output_unit,'(2(A,F8.5))') routinen//' filter(SC): occ_in',occ_in,' occ_out',occ_out 00591 ENDIF 00592 ELSE 00593 CALL cp_dbcsr_copy(SC,C_OLD,error=error) 00594 ENDIF 00595 ! 00596 ! compute P = C'*SC 00597 CALL cp_dbcsr_multiply('T','N',1.0_dp,C_OLD,SC,0.0_dp,P,error=error) 00598 IF (qs_ot_env1%settings%eps_irac_filter_matrix.GT.0.0_dp) THEN 00599 occ_in = cp_dbcsr_get_occupation(p) 00600 CALL cp_dbcsr_filter(p,qs_ot_env1%settings%eps_irac_filter_matrix,error=error) 00601 occ_out = cp_dbcsr_get_occupation(p) 00602 IF(output_unit>0.AND.qs_ot_env1%settings%eps_irac_filter_matrix.GT.0.0_dp) & 00603 WRITE(output_unit,'(2(A,F8.5))') routinen//' filter(P): occ_in',occ_in,' occ_out',occ_out 00604 ENDIF 00605 ! 00606 ! check ||P-1||_f and ||P-1||_gct 00607 CALL cp_dbcsr_add_on_diag(P, -1.0_dp, error=error) 00608 norm_fro = cp_dbcsr_frobenius_norm(P) 00609 norm_gct = cp_dbcsr_gershgorin_norm(P) 00610 CALL cp_dbcsr_add_on_diag(P, 1.0_dp, error=error) 00611 norm = MIN(norm_gct,norm_fro) 00612 CALL qs_ot_ref_decide(qs_ot_env1,norm,ortho_irac) 00613 IF(output_unit>0) WRITE(output_unit,'(A,I3,A,E12.5,A)') & 00614 & routineN,0,': ||P-I||=',norm,& 00615 & ', ortho_irac = '//ortho_irac 00616 ! 00617 ! select the orthogonality method 00618 SELECT CASE(ortho_irac) 00619 CASE("CHOL") 00620 CALL qs_ot_ref_chol(qs_ot_env,C_OLD,C_TMP,C_NEW,P,SC,update,error) 00621 CASE("LWDN") 00622 CALL qs_ot_ref_lwdn(qs_ot_env,C_OLD,C_TMP,C_NEW,P,SC,update,error) 00623 CASE("POLY") 00624 CALL qs_ot_ref_poly(qs_ot_env,C_OLD,C_TMP,C_NEW,P,SC,norm,update,output_unit,error) 00625 CASE DEFAULT 00626 CALL stop_program(routineN,moduleN,__LINE__,"Wrong argument") 00627 END SELECT 00628 ! 00629 ! We update the C_i+1 and localization 00630 IF(update) THEN 00631 IF(on_the_fly_loc) THEN 00632 IF(output_unit>0) WRITE(output_unit,'(A)') & 00633 & routineN//' we localize C' 00634 CALL qs_ot_on_the_fly_localize(qs_ot_env,C_NEW,SC,G_OLD,D,error) 00635 ENDIF 00636 CALL cp_dbcsr_copy(C_OLD,C_NEW,error=error) 00637 ENDIF 00638 ! 00639 CALL timestop(handle) 00640 END SUBROUTINE qs_ot_get_orbitals_ref 00641 00642 SUBROUTINE qs_ot_refine(P,FY,P2,T,irac_degree,eps_irac_filter_matrix,output_unit,error) 00643 !---------------------------------------------------------------------- 00644 ! refinement polynomial of degree 2,3 and 4 (PRB 70, 193102 (2004)) 00645 !---------------------------------------------------------------------- 00646 00647 TYPE(cp_dbcsr_type), INTENT(inout) :: P, FY, P2, T 00648 INTEGER, INTENT(in) :: irac_degree 00649 REAL(dp), INTENT(in) :: eps_irac_filter_matrix 00650 INTEGER, INTENT(in) :: output_unit 00651 TYPE(cp_error_type), INTENT(inout) :: error 00652 00653 CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_refine', 00654 routineP = moduleN//':'//routineN 00655 REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp 00656 00657 INTEGER :: handle, k 00658 REAL(dp) :: occ_in, occ_out, r 00659 00660 CALL timeset(routineN,handle) 00661 00662 CALL cp_dbcsr_get_info(P,nfullcols_total=k) 00663 SELECT CASE(irac_degree) 00664 CASE(2) 00665 ! C_out = C_in * ( 15/8 * I - 10/8 * P + 3/8 * P^2) 00666 r = 3.0_dp/8.0_dp 00667 CALL cp_dbcsr_multiply('N', 'N', r, P, P, 0.0_dp, FY, error=error) 00668 IF (eps_irac_filter_matrix.GT.0.0_dp) THEN 00669 occ_in = cp_dbcsr_get_occupation(fy) 00670 CALL cp_dbcsr_filter(fy,eps_irac_filter_matrix,error=error) 00671 occ_out = cp_dbcsr_get_occupation(fy) 00672 IF(output_unit>0.AND.eps_irac_filter_matrix.GT.0.0_dp) & 00673 WRITE(output_unit,'(2(A,F8.5))') routinen//' filter(FY): occ_in',occ_in,' occ_out',occ_out 00674 ENDIF 00675 r = -10.0_dp/8.0_dp 00676 CALL cp_dbcsr_add(FY, P, alpha_scalar=1.0_dp, beta_scalar=r, error=error) 00677 r = 15.0_dp/8.0_dp 00678 CALL cp_dbcsr_add_on_diag(FY, alpha_scalar=r, error=error) 00679 CASE(3) 00680 ! C_out = C_in * ( 35/16 * I - 35/16 * P + 21/16 * P^2 - 5/16 P^3) 00681 CALL cp_dbcsr_multiply('N', 'N', 1.0_dp, P, P, 0.0_dp, P2, error=error) 00682 IF (eps_irac_filter_matrix.GT.0.0_dp) THEN 00683 occ_in = cp_dbcsr_get_occupation(p2) 00684 CALL cp_dbcsr_filter(p2,eps_irac_filter_matrix,error=error) 00685 occ_out = cp_dbcsr_get_occupation(p2) 00686 IF(output_unit>0.AND.eps_irac_filter_matrix.GT.0.0_dp) & 00687 WRITE(output_unit,'(2(A,F8.5))') routinen//' filter(P2): occ_in',occ_in,' occ_out',occ_out 00688 ENDIF 00689 r = -5.0_dp/16.0_dp 00690 CALL cp_dbcsr_multiply('N', 'N', r, P2, P, 0.0_dp, FY, error=error) 00691 IF (eps_irac_filter_matrix.GT.0.0_dp) THEN 00692 occ_in = cp_dbcsr_get_occupation(fy) 00693 CALL cp_dbcsr_filter(fy,eps_irac_filter_matrix,error=error) 00694 occ_out = cp_dbcsr_get_occupation(fy) 00695 IF(output_unit>0.AND.eps_irac_filter_matrix.GT.0.0_dp) & 00696 WRITE(output_unit,'(2(A,F8.5))') routinen//' filter(FY): occ_in',occ_in,' occ_out',occ_out 00697 ENDIF 00698 r = 21.0_dp/16.0_dp 00699 CALL cp_dbcsr_add(FY, P2, alpha_scalar=1.0_dp, beta_scalar=r, error=error) 00700 r = -35.0_dp/16.0_dp 00701 CALL cp_dbcsr_add(FY, P, alpha_scalar=1.0_dp, beta_scalar=r, error=error) 00702 r = 35.0_dp/16.0_dp 00703 CALL cp_dbcsr_add_on_diag(FY, alpha_scalar=r, error=error) 00704 CASE(4) 00705 ! C_out = C_in * ( 315/128 * I - 420/128 * P + 378/128 * P^2 - 180/128 P^3 + 35/128 P^4 ) 00706 ! = C_in * ( 315/128 * I - 420/128 * P + 378/128 * P^2 + ( - 180/128 * P + 35/128 * P^2 ) * P^2 ) 00707 CALL cp_dbcsr_multiply('N', 'N', 1.0_dp, P, P, 0.0_dp, P2, error=error) ! P^2 00708 IF (eps_irac_filter_matrix.GT.0.0_dp) THEN 00709 occ_in = cp_dbcsr_get_occupation(p2) 00710 CALL cp_dbcsr_filter(p2,eps_irac_filter_matrix,error=error) 00711 occ_out = cp_dbcsr_get_occupation(p2) 00712 IF(output_unit>0.AND.eps_irac_filter_matrix.GT.0.0_dp) & 00713 WRITE(output_unit,'(2(A,F8.5))') routinen//' filter(P2): occ_in',occ_in,' occ_out',occ_out 00714 ENDIF 00715 r = -180.0_dp/128.0_dp 00716 CALL cp_dbcsr_add(T, P, alpha_scalar=0.0_dp, beta_scalar=r, error=error) ! T=-180/128*P 00717 r = 35.0_dp/128.0_dp 00718 CALL cp_dbcsr_add(T, P2, alpha_scalar=1.0_dp, beta_scalar=r, error=error) ! T=T+35/128*P^2 00719 CALL cp_dbcsr_multiply('N', 'N', 1.0_dp, T, P2, 0.0_dp, FY, error=error) ! Y=T*P^2 00720 IF (eps_irac_filter_matrix.GT.0.0_dp) THEN 00721 occ_in = cp_dbcsr_get_occupation(fy) 00722 CALL cp_dbcsr_filter(fy,eps_irac_filter_matrix,error=error) 00723 occ_out = cp_dbcsr_get_occupation(fy) 00724 IF(output_unit>0.AND.eps_irac_filter_matrix.GT.0.0_dp) & 00725 WRITE(output_unit,'(2(A,F8.5))') routinen//' filter(FY): occ_in',occ_in,' occ_out',occ_out 00726 ENDIF 00727 r = 378.0_dp/128.0_dp 00728 CALL cp_dbcsr_add(FY, P2, alpha_scalar=1.0_dp, beta_scalar=r, error=error)! Y=Y+378/128*P^2 00729 r = -420.0_dp/128.0_dp 00730 CALL cp_dbcsr_add(FY, P, alpha_scalar=1.0_dp, beta_scalar=r, error=error) ! Y=Y-420/128*P 00731 r = 315.0_dp/128.0_dp 00732 CALL cp_dbcsr_add_on_diag(FY, alpha_scalar=r, error=error) ! Y=Y+315/128*I 00733 CASE DEFAULT 00734 CALL stop_program(routineN,moduleN,__LINE__,"This irac_order NYI") 00735 END SELECT 00736 CALL timestop(handle) 00737 END SUBROUTINE qs_ot_refine 00738 00739 00740 ! ***************************************************************************** 00741 SUBROUTINE qs_ot_get_derivative_ref(matrix_hc,matrix_x,matrix_sx,matrix_gx, & 00742 & qs_ot_env,output_unit,error) 00743 TYPE(cp_dbcsr_type), POINTER :: matrix_hc, matrix_x, 00744 matrix_sx, matrix_gx 00745 TYPE(qs_ot_type) :: qs_ot_env 00746 INTEGER, INTENT(IN) :: output_unit 00747 TYPE(cp_error_type), INTENT(inout) :: error 00748 00749 CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative_ref', 00750 routineP = moduleN//':'//routineN 00751 REAL(KIND=dp), PARAMETER :: rone = 1.0_dp , rzero = 0.0_dp 00752 00753 INTEGER :: handle, k, n 00754 LOGICAL :: mixed_precision 00755 REAL(dp) :: occ_in, occ_out 00756 TYPE(cp_dbcsr_type), POINTER :: C, CHC, G, G_dp, HC, SC 00757 00758 CALL timeset(routineN,handle) 00759 00760 mixed_precision = qs_ot_env%settings%mixed_precision 00761 CALL cp_dbcsr_get_info(matrix_x,nfullrows_total=n,nfullcols_total=k) 00762 ! 00763 C => matrix_x ! NBsf*NOcc 00764 SC => matrix_sx ! NBsf*NOcc need to be up2date 00765 HC => matrix_hc ! NBsf*NOcc 00766 G => matrix_gx ! NBsf*NOcc 00767 CHC => qs_ot_env%buf1_k_k_sym ! buffer 00768 G_dp => qs_ot_env%buf1_n_k_dp ! buffer dp 00769 00770 IF(mixed_precision) THEN 00771 ! C'*(H*C) 00772 CALL cp_dbcsr_multiply('T','N',1.0_dp,C,HC,rzero,CHC,error=error) 00773 IF (qs_ot_env%settings%eps_irac_filter_matrix.GT.0.0_dp) THEN 00774 occ_in = cp_dbcsr_get_occupation(chc) 00775 CALL cp_dbcsr_filter(chc,qs_ot_env%settings%eps_irac_filter_matrix,error=error) 00776 occ_out = cp_dbcsr_get_occupation(chc) 00777 IF(output_unit>0.AND.qs_ot_env%settings%eps_irac_filter_matrix.GT.0.0_dp) & 00778 WRITE(output_unit,'(2(A,F8.5))') routinen//' filter(CHC): occ_in',occ_in,' occ_out',occ_out 00779 ENDIF 00780 ! (S*C)*(C'*H*C) 00781 CALL cp_dbcsr_multiply('N','N',1.0_dp,SC,CHC,0.0_dp,G_dp,error=error) 00782 IF (qs_ot_env%settings%eps_irac_filter_matrix.GT.0.0_dp) THEN 00783 occ_in = cp_dbcsr_get_occupation(g_dp) 00784 CALL cp_dbcsr_filter(g_dp,qs_ot_env%settings%eps_irac_filter_matrix,error=error) 00785 occ_out = cp_dbcsr_get_occupation(g_dp) 00786 IF(output_unit>0.AND.qs_ot_env%settings%eps_irac_filter_matrix.GT.0.0_dp) & 00787 WRITE(output_unit,'(2(A,F8.5))') routinen//' filter(G_dp): occ_in',occ_in,' occ_out',occ_out 00788 ENDIF 00789 ! G = 2*(1-S*C*C')*H*C 00790 CALL cp_dbcsr_add(G_dp,HC,alpha_scalar=-1.0_dp,beta_scalar=1.0_dp,error=error) 00791 CALL cp_dbcsr_copy(G,G_dp,error=error) 00792 ELSE 00793 ! C'*(H*C) 00794 CALL cp_dbcsr_multiply('T','N',1.0_dp,C,HC,0.0_dp,CHC,error=error) 00795 IF (qs_ot_env%settings%eps_irac_filter_matrix.GT.0.0_dp) THEN 00796 occ_in = cp_dbcsr_get_occupation(chc) 00797 CALL cp_dbcsr_filter(chc,qs_ot_env%settings%eps_irac_filter_matrix,error=error) 00798 occ_out = cp_dbcsr_get_occupation(chc) 00799 IF(output_unit>0.AND.qs_ot_env%settings%eps_irac_filter_matrix.GT.0.0_dp) & 00800 WRITE(output_unit,'(2(A,F8.5))') routinen//' filter(CHC): occ_in',occ_in,' occ_out',occ_out 00801 ENDIF 00802 ! (S*C)*(C'*H*C) 00803 CALL cp_dbcsr_multiply('N','N',1.0_dp,SC,CHC,0.0_dp,G,error=error) 00804 IF (qs_ot_env%settings%eps_irac_filter_matrix.GT.0.0_dp) THEN 00805 occ_in = cp_dbcsr_get_occupation(g) 00806 CALL cp_dbcsr_filter(g,qs_ot_env%settings%eps_irac_filter_matrix,error=error) 00807 occ_out = cp_dbcsr_get_occupation(g) 00808 IF(output_unit>0.AND.qs_ot_env%settings%eps_irac_filter_matrix.GT.0.0_dp) & 00809 WRITE(output_unit,'(2(A,F8.5))') routinen//' filter(G): occ_in',occ_in,' occ_out',occ_out 00810 ENDIF 00811 ! G = 2*(1-S*C*C')*H*C 00812 CALL cp_dbcsr_add(G,HC,alpha_scalar=-1.0_dp,beta_scalar=1.0_dp,error=error) 00813 ENDIF 00814 ! 00815 CALL timestop(handle) 00816 END SUBROUTINE qs_ot_get_derivative_ref 00817 00818 ! computes p=x*S*x and the matrix functionals related matrices 00819 ! ***************************************************************************** 00820 SUBROUTINE qs_ot_get_p(matrix_x,matrix_sx,qs_ot_env,error) 00821 00822 TYPE(cp_dbcsr_type), POINTER :: matrix_x, matrix_sx 00823 TYPE(qs_ot_type) :: qs_ot_env 00824 TYPE(cp_error_type), INTENT(inout) :: error 00825 00826 CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_p', 00827 routineP = moduleN//':'//routineN 00828 REAL(KIND=dp), PARAMETER :: rone = 1.0_dp , rzero = 0.0_dp 00829 00830 INTEGER :: handle, k, n 00831 00832 CALL timeset(routineN,handle) 00833 00834 CALL cp_dbcsr_get_info(matrix_x,nfullrows_total=n,nfullcols_total=k) 00835 00836 ! get the overlap 00837 CALL cp_dbcsr_multiply('T','N',rone,matrix_x,matrix_sx,rzero,& 00838 qs_ot_env%matrix_p,error=error) 00839 ! get an upper bound for the largest eigenvalue 00840 qs_ot_env % largest_eval_upper_bound = cp_dbcsr_gershgorin_norm(qs_ot_env%matrix_p) 00841 CALL decide_strategy(qs_ot_env) 00842 IF (qs_ot_env % do_taylor) THEN 00843 CALL qs_ot_p2m_taylor(qs_ot_env,error=error) 00844 ELSE 00845 CALL qs_ot_p2m_diag(qs_ot_env,error=error) 00846 ENDIF 00847 00848 IF (qs_ot_env % settings % do_rotation) THEN 00849 CALL qs_ot_generate_rotation(qs_ot_env,error=error) 00850 ENDIF 00851 00852 CALL timestop(handle) 00853 00854 END SUBROUTINE qs_ot_get_p 00855 00856 ! ***************************************************************************** 00863 SUBROUTINE qs_ot_generate_rotation(qs_ot_env,error) 00864 00865 TYPE(qs_ot_type) :: qs_ot_env 00866 TYPE(cp_error_type), INTENT(inout) :: error 00867 00868 CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_generate_rotation', 00869 routineP = moduleN//':'//routineN 00870 COMPLEX(KIND=dp), PARAMETER :: cone = (1.0_dp,0.0_dp), 00871 czero = (0.0_dp,0.0_dp) 00872 00873 COMPLEX(KIND=dp), ALLOCATABLE, 00874 DIMENSION(:) :: evals_exp 00875 COMPLEX(KIND=dp), DIMENSION(:), POINTER :: data_z 00876 INTEGER :: blk, col, handle, k, row 00877 LOGICAL :: found 00878 REAL(KIND=dp), DIMENSION(:), POINTER :: data_d 00879 TYPE(cp_dbcsr_iterator) :: iter 00880 TYPE(cp_dbcsr_type) :: cmat_u, cmat_x 00881 00882 CALL timeset(routineN,handle) 00883 00884 CALL cp_dbcsr_get_info(qs_ot_env%rot_mat_x,nfullrows_total=k) 00885 00886 IF (k/=0) THEN 00887 CALL cp_dbcsr_init(cmat_x, error=error) 00888 CALL cp_dbcsr_init(cmat_u, error=error) 00889 CALL cp_dbcsr_copy(cmat_x,qs_ot_env%rot_mat_evec,name='cmat_x',error=error) 00890 CALL cp_dbcsr_copy(cmat_u,qs_ot_env%rot_mat_evec,name='cmat_u',error=error) 00891 ALLOCATE(evals_exp(k)) 00892 00893 ! rot_mat_u = exp(rot_mat_x) 00894 ! i rot_mat_x is hermitian, so go over the complex variables for diag 00895 !vwCALL cp_cfm_get_info(cmat_x,local_data=local_data_c,error=error) 00896 !vwCALL cp_fm_get_info(qs_ot_env%rot_mat_x,local_data=local_data_r,error=error) 00897 !vwlocal_data_c=CMPLX(0.0_dp,local_data_r,KIND=dp) 00898 CALL cp_dbcsr_iterator_start(iter, cmat_x) 00899 DO WHILE (cp_dbcsr_iterator_blocks_left(iter)) 00900 CALL cp_dbcsr_iterator_next_block(iter, row, col, data_z, blk) 00901 CALL cp_dbcsr_get_block_p(qs_ot_env%rot_mat_x, row, col, data_d, found) 00902 IF(.NOT.found) THEN 00903 WRITE(*,*) routineN//' .NOT.found' 00904 !stop 00905 ELSE 00906 data_z=CMPLX(0.0_dp,data_d,KIND=dp) 00907 ENDIF 00908 ENDDO 00909 CALL cp_dbcsr_iterator_stop(iter) 00910 00911 00912 CALL cp_dbcsr_heevd(cmat_x,qs_ot_env%rot_mat_evec,qs_ot_env%rot_mat_evals,& 00913 qs_ot_env%para_env, qs_ot_env%blacs_env, error=error) 00914 evals_exp(:)=EXP( (0.0_dp,-1.0_dp) * qs_ot_env%rot_mat_evals(:) ) 00915 CALL cp_dbcsr_copy(cmat_x,qs_ot_env%rot_mat_evec,error=error) 00916 CALL cp_dbcsr_scale_by_vector(cmat_x,alpha=evals_exp,side='right',error=error) 00917 CALL cp_dbcsr_multiply('N','C',cone,cmat_x,qs_ot_env%rot_mat_evec,czero,cmat_u,error=error) 00918 CALL cp_dbcsr_copy(qs_ot_env%rot_mat_u, cmat_u, keep_imaginary=.FALSE., error=error) 00919 CALL cp_dbcsr_release(cmat_x, error=error) 00920 CALL cp_dbcsr_release(cmat_u, error=error) 00921 DEALLOCATE(evals_exp) 00922 END IF 00923 00924 CALL timestop(handle) 00925 00926 END SUBROUTINE qs_ot_generate_rotation 00927 00928 ! ***************************************************************************** 00935 SUBROUTINE qs_ot_rot_mat_derivative(qs_ot_env,error) 00936 TYPE(qs_ot_type) :: qs_ot_env 00937 TYPE(cp_error_type), INTENT(inout) :: error 00938 00939 CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_rot_mat_derivative', 00940 routineP = moduleN//':'//routineN 00941 00942 COMPLEX(KIND=dp), PARAMETER :: cI = (0.0_dp,1.0_dp), 00943 cone = (1.0_dp,0.0_dp), 00944 czero = (0.0_dp,0.0_dp) 00945 00946 INTEGER :: handle, i, j, k 00947 REAL(KIND=dp) :: e1, e2 00948 REAL(KIND=dp), DIMENSION(:, :), POINTER :: data_d 00949 TYPE(cp_dbcsr_type) :: cmat_buf1, cmat_buf2 00950 TYPE(cp_dbcsr_iterator) :: iter 00951 COMPLEX(dp), DIMENSION(:,:), POINTER :: data_z 00952 INTEGER::row, col, blk,row_offset, col_offset,row_size, col_size 00953 LOGICAL::found 00954 00955 CALL timeset(routineN,handle) 00956 00957 CALL cp_dbcsr_get_info(qs_ot_env%rot_mat_u,nfullrows_total=k) 00958 IF (k/=0) THEN 00959 CALL cp_dbcsr_copy(qs_ot_env%matrix_buf1,qs_ot_env%rot_mat_dedu,error=error) 00960 ! now we get to the derivative wrt the antisymmetric matrix rot_mat_x 00961 CALL cp_dbcsr_init(cmat_buf1, error=error) 00962 CALL cp_dbcsr_init(cmat_buf2, error=error) 00963 CALL cp_dbcsr_copy(cmat_buf1,qs_ot_env%rot_mat_evec,"cmat_buf1",error=error) 00964 CALL cp_dbcsr_copy(cmat_buf2,qs_ot_env%rot_mat_evec,"cmat_buf2",error=error) 00965 00966 ! init cmat_buf1 00967 !CALL cp_fm_get_info(qs_ot_env%matrix_buf1,matrix_struct=fm_struct, local_data=local_data_r,error=error) 00968 !CALL cp_cfm_get_info(cmat_buf1, nrow_local=nrow_local, ncol_local=ncol_local, & 00969 ! row_indices=row_indices, col_indices=col_indices, & 00970 ! local_data=local_data_c,error=error) 00971 !local_data_c=local_data_r 00972 00973 CALL cp_dbcsr_iterator_start(iter, cmat_buf1) 00974 DO WHILE (cp_dbcsr_iterator_blocks_left(iter)) 00975 CALL cp_dbcsr_iterator_next_block(iter, row, col, data_z, blk) 00976 CALL cp_dbcsr_get_block_p(qs_ot_env%matrix_buf1, row, col, data_d, found) 00977 data_z=data_d 00978 ENDDO 00979 CALL cp_dbcsr_iterator_stop(iter) 00980 00981 CALL cp_dbcsr_multiply('T','N',cone,cmat_buf1,qs_ot_env%rot_mat_evec,& 00982 czero,cmat_buf2,error=error) 00983 CALL cp_dbcsr_multiply('C','N',cone,qs_ot_env%rot_mat_evec,cmat_buf2,& 00984 czero,cmat_buf1,error=error) 00985 00986 CALL cp_dbcsr_iterator_start(iter, cmat_buf1) 00987 DO WHILE (cp_dbcsr_iterator_blocks_left (iter)) 00988 CALL cp_dbcsr_iterator_next_block(iter, row, col, data_z,blk,& 00989 row_size=row_size, col_size=col_size, & 00990 row_offset=row_offset, col_offset=col_offset) 00991 DO j=1,col_size 00992 DO i=1,row_size 00993 e1=qs_ot_env%rot_mat_evals(row_offset+i-1) 00994 e2=qs_ot_env%rot_mat_evals(col_offset+j-1) 00995 data_z(i,j)=data_z(i,j)*cint(e1,e2) 00996 ENDDO 00997 ENDDO 00998 ENDDO 00999 CALL cp_dbcsr_iterator_stop(iter) 01000 01001 CALL cp_dbcsr_multiply('N','N',cone,qs_ot_env%rot_mat_evec,cmat_buf1,& 01002 czero,cmat_buf2,error=error) 01003 CALL cp_dbcsr_multiply('N','C',cone,cmat_buf2,qs_ot_env%rot_mat_evec,& 01004 czero,cmat_buf1,error=error) 01005 01006 CALL cp_dbcsr_copy(qs_ot_env%matrix_buf1,cmat_buf1,error=error) 01007 01008 CALL cp_dbcsr_transposed(qs_ot_env%matrix_buf2,qs_ot_env%matrix_buf1,& 01009 shallow_data_copy=.FALSE., use_distribution=cp_dbcsr_distribution(qs_ot_env%matrix_buf3), & 01010 transpose_distribution=.FALSE.,error=error) 01011 CALL cp_dbcsr_add(qs_ot_env%matrix_buf1,qs_ot_env%matrix_buf2,& 01012 alpha_scalar=-1.0_dp,beta_scalar=+1.0_dp,error=error) 01013 CALL cp_dbcsr_copy(qs_ot_env%rot_mat_gx,qs_ot_env%matrix_buf1,error=error) 01014 CALL cp_dbcsr_release(cmat_buf1, error=error) 01015 CALL cp_dbcsr_release(cmat_buf2, error=error) 01016 END IF 01017 CALL timestop(handle) 01018 CONTAINS 01019 ! ***************************************************************************** 01020 FUNCTION cint(e1,e2) 01021 REAL(KIND=dp) :: e1, e2 01022 COMPLEX(KIND=dp) :: cint 01023 01024 COMPLEX(KIND=dp) :: l1, l2, x 01025 INTEGER :: I 01026 01027 l1=(0.0_dp,-1.0_dp)*e1 01028 l2=(0.0_dp,-1.0_dp)*e2 01029 IF (ABS(l1-l2) .GT. 0.5_dp) THEN 01030 cint=(EXP(l1)-EXP(l2))/(l1-l2) 01031 ELSE 01032 x=1.0_dp 01033 cint=0.0_dp 01034 DO I=1,16 01035 cint=cint+x 01036 x=x*(l1-l2)/REAL(I+1,KIND=dp) 01037 ENDDO 01038 cint=cint*EXP(l2) 01039 ENDIF 01040 END FUNCTION cint 01041 END SUBROUTINE qs_ot_rot_mat_derivative 01042 01043 ! 01044 ! decide strategy 01045 ! tries to decide if the taylor expansion of cos(sqrt(xsx)) converges rapidly enough 01046 ! to make a taylor expansion of the functions cos(sqrt(xsx)) and sin(sqrt(xsx))/sqrt(xsx) 01047 ! and their derivatives faster than their computation based on diagonalization 01048 ! since xsx can be very small, especially during dynamics, only a few terms might indeed be needed 01049 ! we find the necessary order N to have largest_eval_upper_bound**(N+1)/(2(N+1))! < eps_taylor 01050 ! 01051 ! ***************************************************************************** 01052 SUBROUTINE decide_strategy(qs_ot_env) 01053 TYPE(qs_ot_type) :: qs_ot_env 01054 01055 INTEGER :: N 01056 REAL(KIND=dp) :: num_error 01057 01058 qs_ot_env % do_taylor = .FALSE. 01059 N=0 01060 num_error=qs_ot_env % largest_eval_upper_bound / ( 2.0_dp ) 01061 DO WHILE (num_error > qs_ot_env % settings % eps_taylor .AND. N <= qs_ot_env % settings % max_taylor) 01062 N=N+1 01063 num_error=num_error * qs_ot_env % largest_eval_upper_bound / REAL(( 2*N+1 )*(2*N+2),KIND=dp) 01064 END DO 01065 qs_ot_env % taylor_order = N 01066 IF ( qs_ot_env % taylor_order <= qs_ot_env % settings % max_taylor) THEN 01067 qs_ot_env % do_taylor = .TRUE. 01068 ENDIF 01069 01070 END SUBROUTINE decide_strategy 01071 01072 ! c=(c0*cos(p^0.5)+x*sin(p^0.5)*p^(-0.5)) x rot_mat_u 01073 ! this assumes that x is already ortho to S*C0, and that p is x*S*x 01074 ! rot_mat_u is an optional rotation matrix 01075 ! ***************************************************************************** 01076 SUBROUTINE qs_ot_get_orbitals(matrix_c, matrix_x, qs_ot_env, error) 01077 01078 TYPE(cp_dbcsr_type), POINTER :: matrix_c, matrix_x 01079 TYPE(qs_ot_type) :: qs_ot_env 01080 TYPE(cp_error_type), INTENT(inout) :: error 01081 01082 CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_orbitals', 01083 routineP = moduleN//':'//routineN 01084 REAL(KIND=dp), PARAMETER :: rone = 1.0_dp , rzero = 0.0_dp 01085 01086 INTEGER :: handle, k, n 01087 TYPE(cp_dbcsr_type), POINTER :: matrix_kk 01088 01089 CALL timeset(routineN,handle) 01090 01091 CALL cp_dbcsr_get_info(matrix_x,nfullrows_total=n,nfullcols_total=k) 01092 01093 ! rotate the multiplying matrices cosp and sinp instead of the result, 01094 ! this should be cheaper for large basis sets 01095 IF (qs_ot_env%settings%do_rotation) THEN 01096 matrix_kk => qs_ot_env%matrix_buf1 01097 CALL cp_dbcsr_multiply('N','N',rone,qs_ot_env%matrix_cosp, & 01098 qs_ot_env%rot_mat_u,rzero,matrix_kk,error=error) 01099 ELSE 01100 matrix_kk => qs_ot_env%matrix_cosp 01101 ENDIF 01102 01103 CALL cp_dbcsr_multiply('N','N',rone,qs_ot_env%matrix_c0,matrix_kk, & 01104 rzero,matrix_c,error=error) 01105 01106 IF (qs_ot_env%settings%do_rotation) THEN 01107 matrix_kk => qs_ot_env%matrix_buf1 01108 CALL cp_dbcsr_multiply('N','N',rone,qs_ot_env%matrix_sinp, & 01109 qs_ot_env%rot_mat_u,rzero,matrix_kk,error=error) 01110 ELSE 01111 matrix_kk => qs_ot_env%matrix_sinp 01112 ENDIF 01113 CALL cp_dbcsr_multiply('N','N',rone,matrix_x,matrix_kk, & 01114 rone ,matrix_c,error=error) 01115 01116 CALL timestop(handle) 01117 01118 END SUBROUTINE qs_ot_get_orbitals 01119 01120 ! ***************************************************************************** 01121 SUBROUTINE qs_ot_get_scp_nddo_coeffs ( qs_ot_env, pscp, error ) 01122 TYPE(qs_ot_type) :: qs_ot_env 01123 TYPE(cp_dbcsr_type), POINTER :: pscp 01124 TYPE(cp_error_type), INTENT(inout) :: error 01125 01126 CALL cp_dbcsr_copy(pscp,qs_ot_env%xmat,error=error) 01127 01128 END SUBROUTINE qs_ot_get_scp_nddo_coeffs 01129 ! ***************************************************************************** 01130 SUBROUTINE qs_ot_get_scp_nddo_derivative ( qs_ot_env, pscp, fscp, error ) 01131 TYPE(qs_ot_type) :: qs_ot_env 01132 TYPE(cp_dbcsr_type), POINTER :: pscp, fscp 01133 TYPE(cp_error_type), INTENT(inout) :: error 01134 01135 CALL cp_dbcsr_copy(qs_ot_env%xmat,pscp,error=error) 01136 CALL cp_dbcsr_copy(qs_ot_env%gxmat,fscp,error=error) 01137 01138 END SUBROUTINE qs_ot_get_scp_nddo_derivative 01139 ! ***************************************************************************** 01140 SUBROUTINE qs_ot_get_scp_dft_coeffs ( qs_ot_env, aux_coeff_set, error ) 01141 TYPE(qs_ot_type) :: qs_ot_env 01142 TYPE(aux_coeff_set_type), POINTER :: aux_coeff_set 01143 TYPE(cp_error_type), INTENT(inout) :: error 01144 01145 INTEGER :: i, icoef, icoef_atom, ikind, 01146 n_els, ncoef_atom, nkind 01147 REAL(dp), DIMENSION(:, :), POINTER :: c 01148 TYPE(aux_coeff_type), POINTER :: local_coeffs 01149 01150 icoef = 0 01151 nkind = SIZE ( aux_coeff_set % coeffs_of_kind ) 01152 DO ikind = 1, nkind 01153 local_coeffs => aux_coeff_set % coeffs_of_kind ( ikind ) % coeffs 01154 IF ( ASSOCIATED ( local_coeffs ) ) THEN 01155 CALL get_aux_coeff ( coeffs = local_coeffs, c = c, & 01156 n_els = n_els, ncoef_atom = ncoef_atom, & 01157 error = error ) 01158 DO i = 1, n_els 01159 DO icoef_atom = 1, ncoef_atom 01160 icoef = icoef + 1 01161 !DBG 01162 ! IF ( icoef == 1 ) & 01163 ! c ( i, icoef_atom ) = qs_ot_env % x ( icoef ) + .05 01164 !DBG 01165 c ( i, icoef_atom ) = qs_ot_env % x ( icoef ) 01166 END DO 01167 END DO 01168 END IF 01169 END DO 01170 END SUBROUTINE qs_ot_get_scp_dft_coeffs 01171 ! this routines sets the SCP derivative to the appropriate 01172 ! qs_ot_env subtype 01173 ! ***************************************************************************** 01174 SUBROUTINE qs_ot_get_scp_dft_derivative ( qs_ot_env, aux_coeff_set, error ) 01175 TYPE(qs_ot_type) :: qs_ot_env 01176 TYPE(aux_coeff_set_type), POINTER :: aux_coeff_set 01177 TYPE(cp_error_type), INTENT(inout) :: error 01178 01179 INTEGER :: i, icoef, icoef_atom, ikind, 01180 n_els, ncoef_atom, nkind 01181 REAL(dp), DIMENSION(:, :), POINTER :: c, fc 01182 TYPE(aux_coeff_type), POINTER :: local_coeffs 01183 01184 icoef = 0 01185 nkind = SIZE ( aux_coeff_set % coeffs_of_kind ) 01186 DO ikind = 1, nkind 01187 local_coeffs => aux_coeff_set % coeffs_of_kind ( ikind ) % coeffs 01188 IF ( ASSOCIATED ( local_coeffs ) ) THEN 01189 CALL get_aux_coeff ( coeffs = local_coeffs, c = c, fc = fc, & 01190 n_els = n_els, ncoef_atom = ncoef_atom, & 01191 error = error ) 01192 DO i = 1, n_els 01193 DO icoef_atom = 1, ncoef_atom 01194 icoef = icoef + 1 01195 qs_ot_env % x ( icoef ) = c ( i, icoef_atom ) 01196 qs_ot_env % gx ( icoef ) = -fc ( i, icoef_atom ) 01197 END DO 01198 END DO 01199 END IF 01200 END DO 01201 END SUBROUTINE qs_ot_get_scp_dft_derivative 01202 01203 ! this routines computes dE/dx=dx, with dx ortho to sc0 01204 ! needs dE/dC=hc,C0,X,SX,p 01205 ! if preconditioned it will not be the derivative, but the lagrangian multiplier 01206 ! is changed so that P*dE/dx is the right derivative (i.e. in the allowed subspace) 01207 ! ***************************************************************************** 01208 SUBROUTINE qs_ot_get_derivative(matrix_hc,matrix_x,matrix_sx,matrix_gx, & 01209 qs_ot_env,error) 01210 TYPE(cp_dbcsr_type), POINTER :: matrix_hc, matrix_x, 01211 matrix_sx, matrix_gx 01212 TYPE(qs_ot_type) :: qs_ot_env 01213 TYPE(cp_error_type), INTENT(inout) :: error 01214 01215 CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative', 01216 routineP = moduleN//':'//routineN 01217 REAL(KIND=dp), PARAMETER :: rone = 1.0_dp , rzero = 0.0_dp 01218 01219 INTEGER :: handle, k, n, ortho_k 01220 TYPE(cp_dbcsr_type), POINTER :: matrix_hc_local, matrix_target 01221 01222 CALL timeset(routineN,handle) 01223 01224 NULLIFY(matrix_hc_local) 01225 01226 CALL cp_dbcsr_get_info(matrix_x,nfullrows_total=n,nfullcols_total=k) 01227 01228 ! could in principle be taken inside qs_ot_get_derivative_* for increased efficiency 01229 ! create a local rotated version of matrix_hc leaving matrix_hc untouched (needed 01230 ! for lagrangian multipliers) 01231 IF (qs_ot_env % settings % do_rotation) THEN 01232 CALL cp_dbcsr_copy(matrix_gx,matrix_hc,error=error) ! use gx as temporary 01233 CALL cp_dbcsr_init_p(matrix_hc_local, error=error) 01234 CALL cp_dbcsr_copy(matrix_hc_local,matrix_hc,name='matrix_hc_local',error=error) 01235 CALL cp_dbcsr_set(matrix_hc_local,0.0_dp,error=error) 01236 CALL cp_dbcsr_multiply('N','T',rone,matrix_gx,qs_ot_env%rot_mat_u,rzero,matrix_hc_local,error=error) 01237 ELSE 01238 matrix_hc_local=>matrix_hc 01239 ENDIF 01240 01241 IF (qs_ot_env % do_taylor) THEN 01242 CALL qs_ot_get_derivative_taylor(matrix_hc_local,matrix_x,matrix_sx,matrix_gx,qs_ot_env,error=error) 01243 ELSE 01244 CALL qs_ot_get_derivative_diag(matrix_hc_local,matrix_x,matrix_sx,matrix_gx,qs_ot_env,error=error) 01245 ENDIF 01246 01247 ! and make it orthogonal 01248 CALL cp_dbcsr_get_info(qs_ot_env%matrix_sc0,nfullcols_total=ortho_k) 01249 01250 IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN 01251 matrix_target => qs_ot_env%matrix_psc0 01252 ELSE 01253 matrix_target => qs_ot_env%matrix_sc0 01254 ENDIF 01255 ! first make the matrix os if not yet valid 01256 IF (.NOT. qs_ot_env%os_valid) THEN 01257 ! this assumes that the preconditioner is a single matrix 01258 ! that maps sc0 onto psc0 01259 01260 IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN 01261 CALL apply_preconditioner(qs_ot_env%preconditioner, qs_ot_env%matrix_sc0, & 01262 qs_ot_env%matrix_psc0 ,error=error) 01263 ENDIF 01264 CALL cp_dbcsr_multiply('T','N',rone,& 01265 qs_ot_env%matrix_sc0,matrix_target, & 01266 rzero,qs_ot_env%matrix_os,& 01267 error=error) 01268 CALL cp_dbcsr_cholesky_decompose(qs_ot_env%matrix_os,& 01269 para_env=qs_ot_env%para_env,blacs_env=qs_ot_env%blacs_env,error=error) 01270 CALL cp_dbcsr_cholesky_invert(qs_ot_env%matrix_os,& 01271 para_env=qs_ot_env%para_env,blacs_env=qs_ot_env%blacs_env,& 01272 upper_to_full=.TRUE.,error=error) 01273 qs_ot_env%os_valid=.TRUE. 01274 ENDIF 01275 CALL cp_dbcsr_multiply('T','N',rone,matrix_target,matrix_gx, & 01276 rzero,qs_ot_env%matrix_buf1_ortho, error=error) 01277 CALL cp_dbcsr_multiply('N','N',rone,qs_ot_env%matrix_os,& 01278 qs_ot_env%matrix_buf1_ortho, rzero,qs_ot_env%matrix_buf2_ortho,error=error) 01279 CALL cp_dbcsr_multiply('N','N',-rone,qs_ot_env%matrix_sc0, & 01280 qs_ot_env%matrix_buf2_ortho, rone,matrix_gx,error=error) 01281 ! also treat the rot_mat gradient here 01282 IF (qs_ot_env%settings%do_rotation) THEN 01283 CALL qs_ot_rot_mat_derivative(qs_ot_env,error=error) 01284 ENDIF 01285 01286 IF (qs_ot_env % settings % do_rotation) THEN 01287 CALL cp_dbcsr_release_p(matrix_hc_local, error=error) 01288 ENDIF 01289 01290 CALL timestop(handle) 01291 01292 END SUBROUTINE qs_ot_get_derivative 01293 01294 ! ***************************************************************************** 01295 SUBROUTINE qs_ot_get_derivative_diag(matrix_hc,matrix_x,matrix_sx,matrix_gx, & 01296 qs_ot_env,error) 01297 01298 TYPE(cp_dbcsr_type), POINTER :: matrix_hc, matrix_x, 01299 matrix_sx, matrix_gx 01300 TYPE(qs_ot_type) :: qs_ot_env 01301 TYPE(cp_error_type), INTENT(inout) :: error 01302 01303 CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative_diag', 01304 routineP = moduleN//':'//routineN 01305 REAL(KIND=dp), PARAMETER :: rone = 1.0_dp , rzero = 0.0_dp 01306 01307 INTEGER :: handle, k, n 01308 01309 CALL timeset(routineN,handle) 01310 01311 CALL cp_dbcsr_get_info(matrix_x,nfullrows_total=n,nfullcols_total=k) 01312 01313 ! go for the derivative now 01314 ! this de/dc*(dX/dx)*sinp 01315 CALL cp_dbcsr_multiply('N','N',rone,matrix_hc,qs_ot_env%matrix_sinp,rzero,matrix_gx,& 01316 error=error) 01317 ! overlap hc*x 01318 CALL cp_dbcsr_multiply('T','N',rone,matrix_hc,matrix_x,rzero,qs_ot_env%matrix_buf2,& 01319 error=error) 01320 ! get it in the basis of the eigenvectors 01321 CALL cp_dbcsr_multiply('N','N',rone,qs_ot_env%matrix_buf2,qs_ot_env%matrix_r,& 01322 rzero,qs_ot_env%matrix_buf1,error=error) 01323 CALL cp_dbcsr_multiply('T','N',rone,qs_ot_env%matrix_r,qs_ot_env%matrix_buf1, & 01324 rzero,qs_ot_env%matrix_buf2,error=error) 01325 01326 ! get the schur product of O_uv*B_uv 01327 CALL cp_dbcsr_hadamard_product(qs_ot_env%matrix_buf2,qs_ot_env%matrix_sinp_b, & 01328 qs_ot_env%matrix_buf3,error=error) 01329 01330 ! overlap hc*c0 01331 CALL cp_dbcsr_multiply('T','N',rone,matrix_hc,qs_ot_env%matrix_c0,rzero, & 01332 qs_ot_env%matrix_buf2,error=error) 01333 ! get it in the basis of the eigenvectors 01334 CALL cp_dbcsr_multiply('N','N',rone,qs_ot_env%matrix_buf2,qs_ot_env%matrix_r, & 01335 rzero,qs_ot_env%matrix_buf1,error=error) 01336 CALL cp_dbcsr_multiply('T','N',rone,qs_ot_env%matrix_r, qs_ot_env%matrix_buf1, & 01337 rzero,qs_ot_env%matrix_buf2,error=error) 01338 01339 CALL cp_dbcsr_hadamard_product(qs_ot_env%matrix_buf2,qs_ot_env%matrix_cosp_b, & 01340 qs_ot_env%matrix_buf4,error=error) 01341 01342 ! add the two bs and compute b+b^T 01343 CALL cp_dbcsr_add(qs_ot_env%matrix_buf3,qs_ot_env%matrix_buf4,& 01344 alpha_scalar=1.0_dp,beta_scalar=1.0_dp,error=error) 01345 01346 ! get the b in the eigenvector basis 01347 CALL cp_dbcsr_multiply('N','T',rone,qs_ot_env%matrix_buf3,qs_ot_env%matrix_r, & 01348 rzero,qs_ot_env%matrix_buf1,error=error) 01349 CALL cp_dbcsr_multiply('N','N',rone,qs_ot_env%matrix_r,qs_ot_env%matrix_buf1, & 01350 rzero,qs_ot_env%matrix_buf3,error=error) 01351 CALL cp_dbcsr_transposed(qs_ot_env%matrix_buf1,qs_ot_env%matrix_buf3,& 01352 shallow_data_copy=.FALSE., use_distribution=cp_dbcsr_distribution(qs_ot_env%matrix_buf3), & 01353 transpose_distribution=.FALSE.,error=error) 01354 CALL cp_dbcsr_add(qs_ot_env%matrix_buf3,qs_ot_env%matrix_buf1,& 01355 alpha_scalar=1.0_dp,beta_scalar=1.0_dp,error=error) 01356 01357 ! and add to the derivative 01358 CALL cp_dbcsr_multiply('N','N',rone,matrix_sx,qs_ot_env%matrix_buf3, & 01359 rone,matrix_gx,error=error) 01360 CALL timestop(handle) 01361 01362 END SUBROUTINE qs_ot_get_derivative_diag 01363 01364 ! compute the derivative of the taylor expansion below 01365 ! ***************************************************************************** 01366 SUBROUTINE qs_ot_get_derivative_taylor(matrix_hc,matrix_x,matrix_sx,matrix_gx, & 01367 qs_ot_env, error) 01368 01369 TYPE(cp_dbcsr_type), POINTER :: matrix_hc, matrix_x, 01370 matrix_sx, matrix_gx 01371 TYPE(qs_ot_type) :: qs_ot_env 01372 TYPE(cp_error_type), INTENT(inout) :: error 01373 01374 CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_get_derivative_taylor', 01375 routineP = moduleN//':'//routineN 01376 REAL(KIND=dp), PARAMETER :: rone = 1.0_dp , rzero = 0.0_dp 01377 01378 INTEGER :: handle, i, k, n 01379 REAL(KIND=dp) :: cosfactor, sinfactor 01380 TYPE(cp_dbcsr_type), POINTER :: matrix_left, matrix_right 01381 01382 CALL timeset(routineN,handle) 01383 01384 CALL cp_dbcsr_get_info(matrix_x,nfullrows_total=n,nfullcols_total=k) 01385 01386 ! go for the derivative now 01387 ! this de/dc*(dX/dx)*sinp i.e. zeroth order 01388 CALL cp_dbcsr_multiply('N','N',rone,matrix_hc,qs_ot_env%matrix_sinp,rzero,matrix_gx,& 01389 error=error) 01390 01391 IF (qs_ot_env % taylor_order .LE. 0) THEN 01392 CALL timestop(handle) 01393 RETURN 01394 ENDIF 01395 01396 ! we store the matrix that will multiply sx in matrix_r 01397 CALL cp_dbcsr_set(qs_ot_env%matrix_r,rzero,error=error) 01398 01399 ! just better names for matrix_cosp_b and matrix_sinp_b (they are buffer space here) 01400 matrix_left => qs_ot_env%matrix_cosp_b 01401 matrix_right => qs_ot_env%matrix_sinp_b 01402 01403 ! overlap hc*x and add its transpose to matrix_left 01404 CALL cp_dbcsr_multiply('T','N',rone,matrix_hc,matrix_x,rzero,matrix_left,& 01405 error=error) 01406 CALL cp_dbcsr_transposed(qs_ot_env%matrix_buf1,matrix_left,& 01407 shallow_data_copy=.FALSE., use_distribution=cp_dbcsr_distribution(matrix_left), & 01408 transpose_distribution=.FALSE., error=error) 01409 CALL cp_dbcsr_add(matrix_left,qs_ot_env%matrix_buf1,& 01410 alpha_scalar=1.0_dp,beta_scalar=1.0_dp,error=error) 01411 CALL cp_dbcsr_copy(matrix_right,matrix_left,error=error) 01412 01413 ! first order 01414 sinfactor=-1.0_dp/(2.0_dp*3.0_dp) 01415 CALL cp_dbcsr_add(qs_ot_env%matrix_r,matrix_left,& 01416 alpha_scalar=1.0_dp,beta_scalar=sinfactor,error=error) 01417 01418 ! M 01419 ! OM+MO 01420 ! OOM+OMO+MOO 01421 ! ... 01422 DO i=2, qs_ot_env % taylor_order 01423 sinfactor=sinfactor * (-1.0_dp)/REAL(2*i * (2*i+1),KIND=dp) 01424 CALL cp_dbcsr_multiply('N','N',rone,qs_ot_env%matrix_p,matrix_left,rzero,qs_ot_env%matrix_buf1,error=error) 01425 CALL cp_dbcsr_multiply('N','N',rone,matrix_right,qs_ot_env%matrix_p,rzero,matrix_left,error=error) 01426 CALL cp_dbcsr_copy(matrix_right,matrix_left,error=error) 01427 CALL cp_dbcsr_add(matrix_left,qs_ot_env%matrix_buf1,& 01428 1.0_dp,1.0_dp,error=error) 01429 CALL cp_dbcsr_add(qs_ot_env%matrix_r,matrix_left,& 01430 alpha_scalar=1.0_dp,beta_scalar=sinfactor,error=error) 01431 ENDDO 01432 01433 ! overlap hc*c0 and add its transpose to matrix_left 01434 CALL cp_dbcsr_multiply('T','N',rone,matrix_hc,qs_ot_env%matrix_c0,rzero,matrix_left,error=error) 01435 CALL cp_dbcsr_transposed(qs_ot_env%matrix_buf1,matrix_left,& 01436 shallow_data_copy=.FALSE., use_distribution=cp_dbcsr_distribution(matrix_left),& 01437 transpose_distribution=.FALSE., error=error) 01438 CALL cp_dbcsr_add(matrix_left,qs_ot_env%matrix_buf1,1.0_dp,1.0_dp,error=error) 01439 CALL cp_dbcsr_copy(matrix_right,matrix_left,error=error) 01440 01441 ! first order 01442 cosfactor=-1.0_dp/(1.0_dp*2.0_dp) 01443 CALL cp_dbcsr_add(qs_ot_env%matrix_r,matrix_left,& 01444 alpha_scalar=1.0_dp,beta_scalar=cosfactor,error=error) 01445 01446 ! M 01447 ! OM+MO 01448 ! OOM+OMO+MOO 01449 ! ... 01450 DO i=2, qs_ot_env % taylor_order 01451 cosfactor=cosfactor * (-1.0_dp)/REAL(2*i * (2*i-1),KIND=dp) 01452 CALL cp_dbcsr_multiply('N','N',rone,qs_ot_env%matrix_p,matrix_left,rzero,qs_ot_env%matrix_buf1,error=error) 01453 CALL cp_dbcsr_multiply('N','N',rone,matrix_right,qs_ot_env%matrix_p,rzero,matrix_left,error=error) 01454 CALL cp_dbcsr_copy(matrix_right,matrix_left,error=error) 01455 CALL cp_dbcsr_add(matrix_left,qs_ot_env%matrix_buf1,1.0_dp,1.0_dp,error=error) 01456 CALL cp_dbcsr_add(qs_ot_env%matrix_r,matrix_left,& 01457 alpha_scalar=1.0_dp,beta_scalar=cosfactor,error=error) 01458 ENDDO 01459 01460 ! and add to the derivative 01461 CALL cp_dbcsr_multiply('N','N',rone,matrix_sx,qs_ot_env%matrix_r,rone,matrix_gx,error=error) 01462 01463 CALL timestop(handle) 01464 01465 END SUBROUTINE qs_ot_get_derivative_taylor 01466 01467 ! computes a taylor expansion. 01468 ! ***************************************************************************** 01469 SUBROUTINE qs_ot_p2m_taylor(qs_ot_env,error) 01470 TYPE(qs_ot_type) :: qs_ot_env 01471 TYPE(cp_error_type), INTENT(inout) :: error 01472 01473 CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_p2m_taylor', 01474 routineP = moduleN//':'//routineN 01475 REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp 01476 01477 INTEGER :: handle, i, k 01478 REAL(KIND=dp) :: cosfactor, sinfactor 01479 01480 CALL timeset(routineN,handle) 01481 01482 ! zeroth order 01483 CALL cp_dbcsr_set(qs_ot_env%matrix_cosp,rzero,error=error) 01484 CALL cp_dbcsr_set(qs_ot_env%matrix_sinp,rzero,error=error) 01485 CALL cp_dbcsr_add_on_diag(qs_ot_env%matrix_cosp,rone,error=error) 01486 CALL cp_dbcsr_add_on_diag(qs_ot_env%matrix_sinp,rone,error=error) 01487 01488 IF (qs_ot_env% taylor_order .LE. 0) THEN 01489 CALL timestop(handle) 01490 RETURN 01491 ENDIF 01492 01493 ! first order 01494 cosfactor=-1.0_dp/(1.0_dp*2.0_dp) 01495 sinfactor=-1.0_dp/(2.0_dp*3.0_dp) 01496 CALL cp_dbcsr_add(qs_ot_env%matrix_cosp,qs_ot_env%matrix_p,alpha_scalar=1.0_dp,beta_scalar=cosfactor,error=error) 01497 CALL cp_dbcsr_add(qs_ot_env%matrix_sinp,qs_ot_env%matrix_p,alpha_scalar=1.0_dp,beta_scalar=sinfactor,error=error) 01498 IF (qs_ot_env% taylor_order .LE. 1) THEN 01499 CALL timestop(handle) 01500 RETURN 01501 ENDIF 01502 01503 ! other orders 01504 CALL cp_dbcsr_get_info(qs_ot_env%matrix_p,nfullrows_total=k) 01505 CALL cp_dbcsr_copy(qs_ot_env%matrix_r,qs_ot_env%matrix_p,error=error) 01506 01507 DO i=2, qs_ot_env%taylor_order 01508 ! new power of p 01509 CALL cp_dbcsr_multiply('N','N',rone,qs_ot_env%matrix_p,qs_ot_env%matrix_r,& 01510 rzero,qs_ot_env%matrix_buf1,error=error) 01511 CALL cp_dbcsr_copy(qs_ot_env%matrix_r,qs_ot_env%matrix_buf1,error=error) 01512 ! add to the taylor expansion so far 01513 cosfactor=cosfactor * (-1.0_dp)/REAL(2*i * (2*i-1),KIND=dp) 01514 sinfactor=sinfactor * (-1.0_dp)/REAL(2*i * (2*i+1),KIND=dp) 01515 CALL cp_dbcsr_add(qs_ot_env%matrix_cosp,qs_ot_env%matrix_r,& 01516 alpha_scalar=1.0_dp,beta_scalar=cosfactor,error=error) 01517 CALL cp_dbcsr_add(qs_ot_env%matrix_sinp,qs_ot_env%matrix_r,& 01518 alpha_scalar=1.0_dp,beta_scalar=sinfactor,error=error) 01519 ENDDO 01520 01521 CALL timestop(handle) 01522 01523 END SUBROUTINE qs_ot_p2m_taylor 01524 01525 ! given p, computes - eigenstuff (matrix_r,evals) 01526 ! - cos(p^0.5),p^(-0.5)*sin(p^0.5) 01527 ! - the real b matrices, needed for the derivatives of these guys 01528 ! cosp_b_ij=(1/(2pii) * int(cos(z^1/2)/((z-eval(i))*(z-eval(j)))) 01529 ! sinp_b_ij=(1/(2pii) * int(z^(-1/2)*sin(z^1/2)/((z-eval(i))*(z-eval(j)))) 01530 ! ***************************************************************************** 01531 SUBROUTINE qs_ot_p2m_diag(qs_ot_env,error) 01532 01533 TYPE(qs_ot_type) :: qs_ot_env 01534 TYPE(cp_error_type), INTENT(inout) :: error 01535 01536 CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_p2m_diag', 01537 routineP = moduleN//':'//routineN 01538 REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp 01539 01540 INTEGER :: blk, col, col_offset, 01541 col_size, handle, i, j, k, 01542 row, row_offset, row_size 01543 REAL(dp), DIMENSION(:, :), POINTER :: DATA 01544 REAL(KIND=dp) :: a, b 01545 TYPE(cp_dbcsr_iterator) :: iter 01546 01547 CALL timeset(routineN,handle) 01548 01549 CALL cp_dbcsr_get_info(qs_ot_env%matrix_p,nfullrows_total=k) 01550 CALL cp_dbcsr_copy(qs_ot_env%matrix_buf1,qs_ot_env%matrix_p,error=error) 01551 CALL cp_dbcsr_syevd(qs_ot_env%matrix_buf1,qs_ot_env%matrix_r,qs_ot_env%evals,& 01552 qs_ot_env%para_env,qs_ot_env%blacs_env,error=error) 01553 DO i=1,k 01554 qs_ot_env%evals(i)=MAX(0.0_dp,qs_ot_env%evals(i)) 01555 ENDDO 01556 01557 !$OMP PARALLEL DO 01558 DO i=1,k 01559 qs_ot_env%dum(i)=COS(SQRT(qs_ot_env%evals(i))) 01560 ENDDO 01561 CALL cp_dbcsr_copy(qs_ot_env%matrix_buf1,qs_ot_env%matrix_r,error=error) 01562 CALL cp_dbcsr_scale_by_vector(qs_ot_env%matrix_buf1,alpha=qs_ot_env%dum,side='right',error=error) 01563 CALL cp_dbcsr_multiply('N','T',rone,qs_ot_env%matrix_r,qs_ot_env%matrix_buf1, & 01564 rzero,qs_ot_env%matrix_cosp,error=error) 01565 01566 !$OMP PARALLEL DO 01567 DO i=1,k 01568 qs_ot_env%dum(i)=qs_ot_sinc(SQRT(qs_ot_env%evals(i))) 01569 ENDDO 01570 CALL cp_dbcsr_copy(qs_ot_env%matrix_buf1,qs_ot_env%matrix_r,error=error) 01571 CALL cp_dbcsr_scale_by_vector(qs_ot_env%matrix_buf1,alpha=qs_ot_env%dum,side='right',error=error) 01572 CALL cp_dbcsr_multiply('N','T',rone,qs_ot_env%matrix_r,qs_ot_env%matrix_buf1, & 01573 rzero,qs_ot_env%matrix_sinp,error=error) 01574 01575 !!$OMP PARALLEL DO PRIVATE(i,j,a,b) 01576 !DO j=1,ncol_local 01577 ! DO i=1,nrow_local 01578 ! a=(SQRT(qs_ot_env%evals(row_indices(i))) & 01579 ! -SQRT(qs_ot_env%evals(col_indices(j))))/2.0_dp 01580 ! b=(SQRT(qs_ot_env%evals(row_indices(i))) & 01581 ! +SQRT(qs_ot_env%evals(col_indices(j))))/2.0_dp 01582 ! qs_ot_env%matrix_cosp_b%local_data(i,j) = -0.5_dp*qs_ot_sinc(a)*qs_ot_sinc(b) 01583 ! ENDDO 01584 !ENDDO 01585 01586 CALL cp_dbcsr_copy(qs_ot_env%matrix_cosp_b,qs_ot_env%matrix_cosp,error=error) 01587 CALL cp_dbcsr_iterator_start(iter, qs_ot_env%matrix_cosp_b) 01588 DO WHILE (cp_dbcsr_iterator_blocks_left (iter)) 01589 CALL cp_dbcsr_iterator_next_block(iter, row, col, DATA,& 01590 block_number=blk, row_size=row_size, col_size=col_size, & 01591 row_offset=row_offset, col_offset=col_offset) 01592 DO j=1,col_size 01593 DO i=1,row_size 01594 a=(SQRT(qs_ot_env%evals( row_offset + i - 1 )) & 01595 -SQRT(qs_ot_env%evals( col_offset + j - 1 )))/2.0_dp 01596 b=(SQRT(qs_ot_env%evals( row_offset + i - 1 )) & 01597 +SQRT(qs_ot_env%evals( col_offset + j - 1 )))/2.0_dp 01598 DATA(i,j) = -0.5_dp*qs_ot_sinc(a)*qs_ot_sinc(b) 01599 ENDDO 01600 ENDDO 01601 ENDDO 01602 CALL cp_dbcsr_iterator_stop(iter) 01603 01604 !!$OMP PARALLEL DO PRIVATE(i,j,a,b) 01605 !DO j=1,ncol_local 01606 ! DO i=1,nrow_local 01607 ! a=SQRT(qs_ot_env%evals(row_indices(i))) 01608 ! b=SQRT(qs_ot_env%evals(col_indices(j))) 01609 ! qs_ot_env%matrix_sinp_b%local_data(i,j)=qs_ot_sincf(a,b) 01610 ! ENDDO 01611 !ENDDO 01612 01613 CALL cp_dbcsr_copy(qs_ot_env%matrix_sinp_b,qs_ot_env%matrix_sinp,error=error) 01614 CALL cp_dbcsr_iterator_start(iter, qs_ot_env%matrix_sinp_b) 01615 DO WHILE (cp_dbcsr_iterator_blocks_left (iter)) 01616 CALL cp_dbcsr_iterator_next_block(iter, row, col, DATA,& 01617 block_number=blk, row_size=row_size, col_size=col_size, & 01618 row_offset=row_offset, col_offset=col_offset) 01619 DO j=1,col_size 01620 DO i=1,row_size 01621 a=SQRT(qs_ot_env%evals( row_offset + i - 1 )) 01622 b=SQRT(qs_ot_env%evals( col_offset + j - 1 )) 01623 DATA(i,j)=qs_ot_sincf(a,b) 01624 ENDDO 01625 ENDDO 01626 ENDDO 01627 CALL cp_dbcsr_iterator_stop(iter) 01628 01629 CALL timestop(handle) 01630 01631 END SUBROUTINE qs_ot_p2m_diag 01632 01633 ! computes sin(x)/x for all values of the argument 01634 ! ***************************************************************************** 01635 FUNCTION qs_ot_sinc(x) 01636 01637 REAL(KIND=dp), INTENT(IN) :: x 01638 REAL(KIND=dp) :: qs_ot_sinc 01639 01640 REAL(KIND=dp), PARAMETER :: q1 = 1.0_dp, q2 = -q1/(2.0_dp *3.0_dp), 01641 q3 = -q2/(4.0_dp *5.0_dp), q4 = -q3/(6.0_dp *7.0_dp), 01642 q5 = -q4/(8.0_dp *9.0_dp), q6 = -q5/(10.0_dp*11.0_dp), 01643 q7 = -q6/(12.0_dp*13.0_dp), q8 = -q7/(14.0_dp*15.0_dp), 01644 q9 = -q8/(16.0_dp*17.0_dp), q10 = -q9/(18.0_dp*19.0_dp) 01645 01646 REAL(KIND=dp) :: y 01647 01648 IF (ABS(x)>0.5_dp) THEN 01649 qs_ot_sinc=SIN(x)/x 01650 ELSE 01651 y=x*x 01652 qs_ot_sinc=q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*(q9+y*(q10))))))))) 01653 ENDIF 01654 END FUNCTION qs_ot_sinc 01655 ! computes (1/(x^2-y^2))*(sinc(x)-sinc(y)) for all positive values of the arguments 01656 ! ***************************************************************************** 01657 FUNCTION qs_ot_sincf(xa,ya) 01658 01659 REAL(KIND=dp), INTENT(IN) :: xa, ya 01660 REAL(KIND=dp) :: qs_ot_sincf 01661 01662 CHARACTER(len=*), PARAMETER :: routineN = 'qs_ot_sincf', 01663 routineP = moduleN//':'//routineN 01664 01665 INTEGER :: i 01666 REAL(KIND=dp) :: a, b, rs, sf, x, xs, y, ybx, 01667 ybxs 01668 01669 ! this is currently a limit of the routine, could be removed rather easily 01670 01671 IF (xa.lt.0) CALL stop_program(routineN,moduleN,__LINE__,"x is negative") 01672 IF (ya.lt.0) CALL stop_program(routineN,moduleN,__LINE__,"y is negative") 01673 01674 IF (xa.lt.ya) THEN 01675 x=ya 01676 y=xa 01677 ELSE 01678 x=xa 01679 y=ya 01680 ENDIF 01681 01682 IF ( x .LT. 0.5_dp ) THEN ! use series, keeping in mind that x,y,x+y,x-y can all be zero 01683 01684 qs_ot_sincf=0.0_dp 01685 IF (x .GT. 0.0_dp) THEN 01686 ybx=y/x 01687 ELSE ! should be irrelevant !? 01688 ybx=0.0_dp 01689 ENDIF 01690 01691 sf=-1.0_dp/((1.0_dp+ybx)*6.0_dp) 01692 rs=1.0_dp 01693 ybxs=ybx 01694 xs=1.0_dp 01695 01696 DO i=1,10 01697 qs_ot_sincf=qs_ot_sincf+sf*rs*xs*(1.0_dp+ybxs) 01698 sf=-sf/(REAL((2*i+2),dp)*REAL((2*i+3),dp)) 01699 rs=rs+ybxs 01700 ybxs=ybxs*ybx 01701 xs=xs*x*x 01702 ENDDO 01703 01704 ELSE ! no series expansion 01705 IF ( x-y .GT. 0.1_dp ) THEN ! safe to use the normal form 01706 qs_ot_sincf=(qs_ot_sinc(x)-qs_ot_sinc(y))/((x+y)*(x-y)) 01707 ELSE 01708 a=(x+y)/2.0_dp 01709 b=(x-y)/2.0_dp ! might be close to zero 01710 ! y (=(a-b)) can not be close to zero since it is close to x>0.5 01711 qs_ot_sincf=(qs_ot_sinc(b)*COS(a)-qs_ot_sinc(a)*COS(b))/(2*x*y) 01712 ENDIF 01713 ENDIF 01714 01715 END FUNCTION qs_ot_sincf 01716 01717 END MODULE qs_ot
1.7.3