|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00012 MODULE qs_ot_minimizer 00013 00014 USE cp_dbcsr_interface, ONLY: cp_dbcsr_add,& 00015 cp_dbcsr_copy,& 00016 cp_dbcsr_get_info,& 00017 cp_dbcsr_scale,& 00018 cp_dbcsr_set,& 00019 cp_dbcsr_trace 00020 USE cp_dbcsr_types, ONLY: cp_dbcsr_p_type,& 00021 cp_dbcsr_type 00022 USE dbcsr_error_handling, ONLY: dbcsr_error_type 00023 USE dbcsr_operations, ONLY: dbcsr_init_random 00024 USE f77_blas 00025 USE kinds, ONLY: dp,& 00026 int_8 00027 USE mathlib, ONLY: diamat_all 00028 USE message_passing, ONLY: mp_sum 00029 USE preconditioner, ONLY: apply_preconditioner 00030 USE qs_ot, ONLY: qs_ot_get_derivative,& 00031 qs_ot_get_derivative_ref,& 00032 qs_ot_get_scp_dft_derivative,& 00033 qs_ot_get_scp_nddo_derivative 00034 USE qs_ot_types, ONLY: qs_ot_type 00035 USE scp_coeff_types, ONLY: aux_coeff_set_type 00036 USE termination, ONLY: stop_program 00037 USE timings, ONLY: timeset,& 00038 timestop 00039 #include "cp_common_uses.h" 00040 00041 IMPLICIT NONE 00042 00043 PRIVATE 00044 00045 PUBLIC :: ot_mini 00046 00047 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ot_minimizer' 00048 00049 CONTAINS 00050 ! 00051 ! the minimizer interface 00052 ! should present all possible modes of minimization 00053 ! these include CG SD DIIS 00054 ! 00055 ! 00056 ! IN the case of nspin != 1 we have a gradient that is distributed over different qs_ot_env. 00057 ! still things remain basically the same, since there are no constraints between the different qs_ot_env 00058 ! we only should take care that the various scalar products are taken over the full vectors. 00059 ! all the information needed and collected can be stored in the fist qs_ot_env only 00060 ! (indicating that the data type for the gradient/position and minization should be separated) 00061 ! 00062 ! ***************************************************************************** 00063 SUBROUTINE ot_mini(qs_ot_env,matrix_hc,output_unit,aux_coeff_set, pscp, fscp, error) 00064 TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env 00065 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00066 POINTER :: matrix_hc 00067 INTEGER, INTENT(IN) :: output_unit 00068 TYPE(aux_coeff_set_type), OPTIONAL, 00069 POINTER :: aux_coeff_set 00070 TYPE(cp_dbcsr_type), OPTIONAL, POINTER :: pscp, fscp 00071 TYPE(cp_error_type), INTENT(inout) :: error 00072 00073 CHARACTER(len=*), PARAMETER :: routineN = 'ot_mini', 00074 routineP = moduleN//':'//routineN 00075 00076 INTEGER :: handle, ispin, nspin 00077 LOGICAL :: do_ener, do_ks, do_scp_dft, 00078 do_scp_nddo 00079 REAL(KIND=dp) :: tmp 00080 00081 CALL timeset(routineN,handle) 00082 00083 nspin=SIZE(qs_ot_env) 00084 00085 do_ks = qs_ot_env ( 1 ) % settings % ks 00086 do_scp_dft = qs_ot_env ( 1 ) % settings % scp_dft 00087 do_scp_nddo = qs_ot_env ( 1 ) % settings % scp_nddo 00088 do_ener = qs_ot_env ( 1 ) % settings % do_ener 00089 00090 qs_ot_env(1)%OT_METHOD_FULL="" 00091 00092 ! compute the gradient for the variables x 00093 IF (.NOT. qs_ot_env(1)%energy_only) THEN 00094 qs_ot_env(1)%gradient=0.0_dp 00095 ! **** SCP 00096 IF ( do_scp_dft ) CALL qs_ot_get_scp_dft_derivative ( qs_ot_env ( 1 ), aux_coeff_set, error ) 00097 IF ( do_scp_nddo ) CALL qs_ot_get_scp_nddo_derivative ( qs_ot_env ( 1 ), pscp, fscp, error ) 00098 ! **** SCP 00099 DO ispin=1,nspin 00100 IF ( do_ks ) THEN 00101 SELECT CASE(qs_ot_env(1)%settings%ot_algorithm) 00102 CASE("TOD") 00103 CALL qs_ot_get_derivative(matrix_hc(ispin)%matrix,qs_ot_env(ispin)%matrix_x, & 00104 qs_ot_env(ispin)%matrix_sx, & 00105 qs_ot_env(ispin)%matrix_gx,qs_ot_env(ispin),error=error) 00106 CASE("REF") 00107 CALL qs_ot_get_derivative_ref(matrix_hc(ispin)%matrix,& 00108 & qs_ot_env(ispin)%matrix_x,qs_ot_env(ispin)%matrix_sx, & 00109 & qs_ot_env(ispin)%matrix_gx,qs_ot_env(ispin),output_unit,error=error) 00110 CASE DEFAULT 00111 CALL stop_program(routineN,moduleN,__LINE__,"ALGORITHM NYI") 00112 END SELECT 00113 END IF 00114 ! and also the gradient along the direction 00115 IF (qs_ot_env(1)%use_dx) THEN 00116 IF ( do_ks ) THEN 00117 CALL cp_dbcsr_trace(qs_ot_env(ispin)%matrix_gx,qs_ot_env(ispin)%matrix_dx,tmp,error=error) 00118 qs_ot_env(1)%gradient=qs_ot_env(1)%gradient+tmp 00119 IF (qs_ot_env(1)%settings%do_rotation) THEN 00120 CALL cp_dbcsr_trace(qs_ot_env(ispin)%rot_mat_gx,qs_ot_env(ispin)%rot_mat_dx,tmp,error=error) 00121 qs_ot_env(1)%gradient=qs_ot_env(1)%gradient+0.5_dp*tmp 00122 ENDIF 00123 END IF 00124 ! ***SCP 00125 IF ( do_scp_dft .AND. ispin == 1 ) THEN 00126 tmp = DOT_PRODUCT ( qs_ot_env ( 1 ) % gx, qs_ot_env ( 1 ) % dx ) 00127 CALL mp_sum ( tmp, qs_ot_env ( 1 ) % scp_para_env % group ) 00128 qs_ot_env ( 1 ) % gradient = qs_ot_env ( 1 ) % gradient + tmp 00129 ENDIF 00130 00131 IF ( do_scp_nddo .AND. ispin == 1 ) THEN 00132 CALL cp_dbcsr_trace ( qs_ot_env ( 1 ) % gxmat, qs_ot_env ( 1 ) % dxmat, tmp, local_sum=.TRUE., & 00133 error=error ) 00134 qs_ot_env ( 1 ) % gradient = qs_ot_env ( 1 ) % gradient + tmp 00135 ENDIF 00136 ! ***SCP 00137 IF (do_ener) THEN 00138 tmp = DOT_PRODUCT ( qs_ot_env ( ispin ) % ener_gx, qs_ot_env ( ispin ) % ener_dx ) 00139 qs_ot_env ( 1 ) % gradient = qs_ot_env ( 1 ) % gradient + tmp 00140 ENDIF 00141 ELSE 00142 IF ( do_ks ) THEN 00143 CALL cp_dbcsr_trace(qs_ot_env(ispin)%matrix_gx,qs_ot_env(ispin)%matrix_gx,tmp,error=error) 00144 qs_ot_env(1)%gradient=qs_ot_env(1)%gradient-tmp 00145 IF (qs_ot_env(1)%settings%do_rotation) THEN 00146 CALL cp_dbcsr_trace(qs_ot_env(ispin)%rot_mat_gx,qs_ot_env(ispin)%rot_mat_gx,tmp,error=error) 00147 qs_ot_env(1)%gradient=qs_ot_env(1)%gradient-0.5_dp*tmp 00148 ENDIF 00149 ENDIF 00150 ! ***SCP 00151 IF ( do_scp_dft .AND. ispin == 1 ) THEN 00152 tmp = DOT_PRODUCT ( qs_ot_env ( 1 ) % gx, qs_ot_env ( 1 ) % gx ) 00153 CALL mp_sum ( tmp, qs_ot_env ( 1 ) % scp_para_env % group ) 00154 qs_ot_env ( 1 ) % gradient = qs_ot_env ( 1 ) % gradient - tmp 00155 ENDIF 00156 00157 IF ( do_scp_nddo .AND. ispin == 1 ) THEN 00158 CALL cp_dbcsr_trace ( qs_ot_env ( 1 ) % gxmat, qs_ot_env ( 1 ) % gxmat, tmp, local_sum=.TRUE.,& 00159 error=error ) 00160 qs_ot_env ( 1 ) % gradient = qs_ot_env ( 1 ) % gradient - tmp 00161 ENDIF 00162 ! ***SCP 00163 IF (do_ener) THEN 00164 tmp = DOT_PRODUCT ( qs_ot_env ( ispin ) % ener_gx, qs_ot_env ( ispin ) % ener_gx ) 00165 qs_ot_env ( 1 ) % gradient = qs_ot_env ( 1 ) % gradient - tmp 00166 ENDIF 00167 ENDIF 00168 ENDDO 00169 ENDIF 00170 00171 SELECT CASE(qs_ot_env(1)%settings%OT_METHOD) 00172 CASE ("CG") 00173 IF (current_point_is_fine(qs_ot_env)) THEN 00174 IF ( ( do_scp_dft .OR. do_scp_nddo ) .AND. .NOT. do_ks ) THEN 00175 qs_ot_env(1)%OT_METHOD_FULL="SCP CG" 00176 ELSEIF ( do_ks .AND. .NOT. do_scp_dft .AND. .NOT. do_scp_nddo ) THEN 00177 qs_ot_env(1)%OT_METHOD_FULL="OT CG" 00178 ELSEIF ( do_ks .AND. ( do_scp_dft .OR. do_scp_nddo ) ) THEN 00179 qs_ot_env(1)%OT_METHOD_FULL="OT CG" 00180 END IF 00181 CALL ot_new_cg_direction(qs_ot_env,error=error) 00182 qs_ot_env(1)%line_search_count=0 00183 ELSE 00184 qs_ot_env(1)%OT_METHOD_FULL="OT LS" 00185 ENDIF 00186 CALL do_line_search(qs_ot_env,error=error) 00187 CASE ("SD") 00188 IF (current_point_is_fine(qs_ot_env)) THEN 00189 IF ( ( do_scp_dft .OR. do_scp_nddo ) .AND. .NOT. do_ks ) THEN 00190 qs_ot_env(1)%OT_METHOD_FULL="SCP SD" 00191 ELSEIF ( do_ks .AND. .NOT. do_scp_dft .AND. .NOT. do_scp_nddo) THEN 00192 qs_ot_env(1)%OT_METHOD_FULL="OT SD" 00193 ELSEIF ( ( do_scp_dft .OR. do_scp_nddo ) .AND. do_ks ) THEN 00194 qs_ot_env(1)%OT_METHOD_FULL="OT SD" 00195 ENDIF 00196 CALL ot_new_sd_direction(qs_ot_env,error=error) 00197 qs_ot_env(1)%line_search_count=0 00198 ELSE 00199 qs_ot_env(1)%OT_METHOD_FULL="OT LS" 00200 ENDIF 00201 CALL do_line_search(qs_ot_env,error=error) 00202 CASE ("DIIS") 00203 IF ( ( do_scp_dft .OR. do_scp_nddo ) .AND. .NOT. do_ks ) THEN 00204 qs_ot_env(1)%OT_METHOD_FULL="SCP DIIS" 00205 ELSEIF ( do_ks .AND. .NOT. do_scp_dft .AND. .NOT. do_scp_nddo ) THEN 00206 qs_ot_env(1)%OT_METHOD_FULL="OT DIIS" 00207 ELSEIF ( ( do_scp_dft .OR. do_scp_nddo ) .AND. do_ks) THEN 00208 qs_ot_env(1)%OT_METHOD_FULL="OT DIIS" 00209 END IF 00210 CALL ot_diis_step(qs_ot_env,error=error) 00211 CASE ("BROY") 00212 qs_ot_env(1)%OT_METHOD_FULL="OT BROY" 00213 CALL ot_broyden_step(qs_ot_env,error=error) 00214 CASE DEFAULT 00215 CALL stop_program(routineN,moduleN,__LINE__,"OT_METHOD NYI") 00216 END SELECT 00217 00218 CALL timestop(handle) 00219 00220 END SUBROUTINE ot_mini 00221 00222 ! 00223 ! checks if the current point is a good point for finding a new direction 00224 ! or if we should improve the line_search, if it is used 00225 ! 00226 ! ***************************************************************************** 00227 FUNCTION current_point_is_fine(qs_ot_env) RESULT(res) 00228 TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env 00229 LOGICAL :: res 00230 00231 res=.FALSE. 00232 00233 ! only if we have a gradient it can be fine 00234 IF (.NOT. qs_ot_env(1)%energy_only ) THEN 00235 00236 ! we have not yet started with the line search 00237 IF (qs_ot_env(1)%line_search_count .EQ. 0) THEN 00238 res=.TRUE. 00239 RETURN 00240 ENDIF 00241 00242 IF (qs_ot_env(1)%line_search_might_be_done) THEN 00243 ! here we put the more complicated logic later 00244 res=.TRUE. 00245 RETURN 00246 ENDIF 00247 00248 ENDIF 00249 00250 END FUNCTION current_point_is_fine 00251 00252 ! 00253 ! performs various kinds of line searches 00254 ! 00255 ! ***************************************************************************** 00256 SUBROUTINE do_line_search(qs_ot_env,error) 00257 TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env 00258 TYPE(cp_error_type), INTENT(inout) :: error 00259 00260 CHARACTER(len=*), PARAMETER :: routineN = 'do_line_search', 00261 routineP = moduleN//':'//routineN 00262 00263 SELECT CASE(qs_ot_env(1)%settings%line_search_method) 00264 CASE("GOLD") 00265 CALL do_line_search_gold(qs_ot_env,error=error) 00266 CASE("3PNT") 00267 CALL do_line_search_3pnt(qs_ot_env,error=error) 00268 CASE("2PNT") 00269 CALL do_line_search_2pnt(qs_ot_env,error=error) 00270 CASE("NONE") 00271 CALL do_line_search_none(qs_ot_env,error=error) 00272 CASE DEFAULT 00273 CALL stop_program(routineN,moduleN,__LINE__,"NYI") 00274 END SELECT 00275 END SUBROUTINE do_line_search 00276 00277 ! ***************************************************************************** 00282 SUBROUTINE take_step(ds,qs_ot_env,error) 00283 REAL(KIND=dp) :: ds 00284 TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env 00285 TYPE(cp_error_type), INTENT(inout) :: error 00286 00287 CHARACTER(len=*), PARAMETER :: routineN = 'take_step', 00288 routineP = moduleN//':'//routineN 00289 00290 INTEGER :: ispin, nspin 00291 LOGICAL :: do_ener, do_ks, do_scp_dft, 00292 do_scp_nddo 00293 00294 nspin=SIZE(qs_ot_env) 00295 00296 do_ks = qs_ot_env ( 1 ) % settings % ks 00297 do_scp_dft = qs_ot_env ( 1 ) % settings % scp_dft 00298 do_scp_nddo = qs_ot_env ( 1 ) % settings % scp_nddo 00299 do_ener = qs_ot_env ( 1 ) % settings % do_ener 00300 00301 ! now update x to take into account this new step 00302 ! either dx or -gx is the direction to use 00303 IF (qs_ot_env(1)%use_dx) THEN 00304 IF ( do_ks ) THEN 00305 DO ispin=1,nspin 00306 CALL cp_dbcsr_add(qs_ot_env(ispin)%matrix_x,qs_ot_env(ispin)%matrix_dx,& 00307 alpha_scalar=1.0_dp,beta_scalar=ds,error=error) 00308 IF (qs_ot_env(ispin)%settings%do_rotation) THEN 00309 CALL cp_dbcsr_add(qs_ot_env(ispin)%rot_mat_x,qs_ot_env(ispin)%rot_mat_dx,& 00310 alpha_scalar=1.0_dp,beta_scalar=ds,error=error) 00311 ENDIF 00312 ENDDO 00313 END IF 00314 ! **** SCP 00315 IF (do_scp_dft) THEN 00316 qs_ot_env(1)%x = qs_ot_env ( 1 ) % x + ds * qs_ot_env ( 1 ) % dx 00317 ENDIF 00318 IF (do_scp_nddo) THEN 00319 CALL cp_dbcsr_add ( qs_ot_env ( 1 ) % xmat, qs_ot_env ( 1 ) % dxmat, & 00320 alpha_scalar=1.0_dp, beta_scalar=ds, error=error ) 00321 ENDIF 00322 ! **** SCP 00323 IF (do_ener) THEN 00324 DO ispin=1,nspin 00325 qs_ot_env(ispin)%ener_x = qs_ot_env ( ispin ) % ener_x + ds * qs_ot_env ( ispin ) % ener_dx 00326 ENDDO 00327 ENDIF 00328 ELSE 00329 IF ( do_ks ) THEN 00330 DO ispin=1,nspin 00331 CALL cp_dbcsr_add(qs_ot_env(ispin)%matrix_x,qs_ot_env(ispin)%matrix_gx,& 00332 alpha_scalar=1.0_dp,beta_scalar=-ds,error=error) 00333 IF (qs_ot_env(ispin)%settings%do_rotation) THEN 00334 CALL cp_dbcsr_add(qs_ot_env(ispin)%rot_mat_x,qs_ot_env(ispin)%rot_mat_gx,& 00335 alpha_scalar=1.0_dp,beta_scalar=-ds,error=error) 00336 ENDIF 00337 ENDDO 00338 ENDIF 00339 ! **** SCP 00340 IF (do_scp_dft) THEN 00341 qs_ot_env(1)%x = qs_ot_env ( 1 ) % x - ds * qs_ot_env ( 1 ) % gx 00342 ENDIF 00343 IF (do_scp_nddo) THEN 00344 CALL cp_dbcsr_add ( qs_ot_env ( 1 ) % xmat, qs_ot_env ( 1 ) % gxmat, & 00345 alpha_scalar=1.0_dp, beta_scalar=-ds, error=error ) 00346 END IF 00347 ! **** SCP 00348 IF (do_ener) THEN 00349 DO ispin=1,nspin 00350 qs_ot_env(ispin)%ener_x = qs_ot_env ( ispin ) % ener_x - ds * qs_ot_env ( ispin ) % ener_gx 00351 ENDDO 00352 ENDIF 00353 ENDIF 00354 END SUBROUTINE take_step 00355 00356 ! implements a golden ratio search as a robust way of minimizing 00357 ! ***************************************************************************** 00358 SUBROUTINE do_line_search_gold(qs_ot_env,error) 00359 00360 TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env 00361 TYPE(cp_error_type), INTENT(inout) :: error 00362 00363 CHARACTER(len=*), PARAMETER :: routineN = 'do_line_search_gold', 00364 routineP = moduleN//':'//routineN 00365 REAL(KIND=dp), PARAMETER :: gold_sec = 0.3819_dp 00366 00367 INTEGER :: count 00368 REAL(KIND=dp) :: ds 00369 00370 ! approx (3-sqrt(5))/2 00371 00372 qs_ot_env(1)%line_search_count=qs_ot_env(1)%line_search_count+1 00373 count=qs_ot_env(1)%line_search_count 00374 qs_ot_env(1)%line_search_might_be_done=.FALSE. 00375 qs_ot_env(1)%energy_only=.TRUE. 00376 00377 IF (count+1 .GT. SIZE(qs_ot_env(1)%OT_pos)) THEN 00378 ! should not happen, we pass with a warning first 00379 ! you can increase the size of OT_pos and the like in qs_ot_env 00380 CALL stop_program(routineN,moduleN,__LINE__,"MAX ITER EXCEEDED : FATAL") 00381 ENDIF 00382 00383 IF (qs_ot_env(1)%line_search_count .EQ. 1) THEN 00384 qs_ot_env(1)%line_search_left = 1 00385 qs_ot_env(1)%line_search_right = 0 00386 qs_ot_env(1)%line_search_mid = 1 00387 qs_ot_env(1)%ot_pos(1) = 0.0_dp 00388 qs_ot_env(1)%ot_energy(1) = qs_ot_env(1)%etotal 00389 qs_ot_env(1)%ot_pos(2) = qs_ot_env(1)%ds_min/gold_sec 00390 ELSE 00391 qs_ot_env(1)%ot_energy(count)=qs_ot_env(1)%etotal 00392 ! it's essentially a book keeping game. 00393 ! keep left on the left, keep (bring) right on the right 00394 ! and mid in between these two 00395 IF (qs_ot_env(1)%line_search_right .EQ. 0) THEN ! we do not yet have the right bracket 00396 IF (qs_ot_env(1)%ot_energy(count-1) .LT. qs_ot_env(1)%ot_energy(count)) THEN 00397 qs_ot_env(1)%line_search_right = count 00398 qs_ot_env(1)%ot_pos(count+1) = qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid)+ & 00399 (qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_right)- & 00400 qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid))*gold_sec 00401 ELSE 00402 qs_ot_env(1)%line_search_left = qs_ot_env(1)%line_search_mid 00403 qs_ot_env(1)%line_search_mid = count 00404 qs_ot_env(1)%ot_pos(count+1) = qs_ot_env(1)%ot_pos(count)/gold_sec ! expand 00405 ENDIF 00406 ELSE 00407 ! first determine where we are and construct the new triplet 00408 IF (qs_ot_env(1)%ot_pos(count) .LT. qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid)) THEN 00409 IF ( qs_ot_env(1)%ot_energy(count) .LT. qs_ot_env(1)%ot_energy(qs_ot_env(1)%line_search_mid))THEN 00410 qs_ot_env(1)%line_search_right = qs_ot_env(1)%line_search_mid 00411 qs_ot_env(1)%line_search_mid = count 00412 ELSE 00413 qs_ot_env(1)%line_search_left = count 00414 ENDIF 00415 ELSE 00416 IF ( qs_ot_env(1)%ot_energy(count) .LT. qs_ot_env(1)%ot_energy(qs_ot_env(1)%line_search_mid))THEN 00417 qs_ot_env(1)%line_search_left = qs_ot_env(1)%line_search_mid 00418 qs_ot_env(1)%line_search_mid = count 00419 ELSE 00420 qs_ot_env(1)%line_search_right = count 00421 ENDIF 00422 ENDIF 00423 ! now find the new point in the largest section 00424 IF ( (qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_right) & 00425 -qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid)) .GT. & 00426 (qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid) & 00427 -qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_left)) ) THEN 00428 qs_ot_env(1)%ot_pos(count+1) = & 00429 qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid) + & 00430 gold_sec*(qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_right) & 00431 -qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid)) 00432 ELSE 00433 qs_ot_env(1)%ot_pos(count+1) = & 00434 qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_left) + & 00435 gold_sec*(qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid) & 00436 -qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_left)) 00437 ENDIF 00438 ! check for termination 00439 IF ( ((qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_right) & 00440 -qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid)) .LT. & 00441 qs_ot_env(1)%ds_min * qs_ot_env(1)%settings%gold_target ) .AND. & 00442 ((qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid) & 00443 -qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_left)).LT. & 00444 qs_ot_env(1)%ds_min * qs_ot_env(1)%settings%gold_target ) ) THEN 00445 qs_ot_env(1)%energy_only=.FALSE. 00446 qs_ot_env(1)%line_search_might_be_done=.TRUE. 00447 ENDIF 00448 ENDIF 00449 ENDIF 00450 ds=qs_ot_env(1)%OT_pos(count+1)-qs_ot_env(1)%OT_pos(count) 00451 qs_ot_env(1)%ds_min=qs_ot_env(1)%OT_pos(count+1) 00452 00453 CALL take_step(ds,qs_ot_env,error=error) 00454 00455 END SUBROUTINE do_line_search_gold 00456 00457 ! ***************************************************************************** 00458 SUBROUTINE do_line_search_3pnt(qs_ot_env,error) 00459 00460 TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env 00461 TYPE(cp_error_type), INTENT(inout) :: error 00462 00463 CHARACTER(len=*), PARAMETER :: routineN = 'do_line_search_3pnt', 00464 routineP = moduleN//':'//routineN 00465 00466 INTEGER :: count 00467 REAL(KIND=dp) :: denom, ds, fa, fb, fc, nom, 00468 pos, val, xa, xb, xc 00469 00470 qs_ot_env(1)%line_search_might_be_done=.FALSE. 00471 qs_ot_env(1)%energy_only=.TRUE. 00472 00473 ! a three point interpolation based on the energy 00474 qs_ot_env(1)%line_search_count=qs_ot_env(1)%line_search_count+1 00475 count=qs_ot_env(1)%line_search_count 00476 qs_ot_env(1)%ot_energy(count)=qs_ot_env(1)%etotal 00477 SELECT CASE(count) 00478 CASE(1) 00479 qs_ot_env(1)%ot_pos(count)=0.0_dp 00480 qs_ot_env(1)%ot_pos(count+1)=qs_ot_env(1)%ds_min*0.8_dp 00481 CASE(2) 00482 IF (qs_ot_env(1)%OT_energy(count).gt.qs_ot_env(1)%OT_energy(count-1)) THEN 00483 qs_ot_env(1)%OT_pos(count+1)=qs_ot_env(1)%ds_min*0.5_dp 00484 ELSE 00485 qs_ot_env(1)%OT_pos(count+1)=qs_ot_env(1)%ds_min*1.4_dp 00486 ENDIF 00487 CASE(3) 00488 xa=qs_ot_env(1)%OT_pos(1) 00489 xb=qs_ot_env(1)%OT_pos(2) 00490 xc=qs_ot_env(1)%OT_pos(3) 00491 fa=qs_ot_env(1)%OT_energy(1) 00492 fb=qs_ot_env(1)%OT_energy(2) 00493 fc=qs_ot_env(1)%OT_energy(3) 00494 nom =(xb-xa)**2*(fb-fc) - (xb-xc)**2*(fb-fa) 00495 denom=(xb-xa)*(fb-fc) - (xb-xc)*(fb-fa) 00496 IF (ABS(denom) .LE. 1.0E-18_dp*MAX(ABS(fb-fc),ABS(fb-fa))) THEN 00497 pos = xb 00498 ELSE 00499 pos = xb-0.5_dp*nom/denom ! position of the stationary point 00500 ENDIF 00501 val = (pos-xa)*(pos-xb)*fc/((xc-xa)*(xc-xb))+ & 00502 (pos-xb)*(pos-xc)*fa/((xa-xb)*(xa-xc))+ & 00503 (pos-xc)*(pos-xa)*fb/((xb-xc)*(xb-xa)) 00504 IF (val.lt.fa .AND. val.le.fb .AND. val.le.fc) THEN ! OK, we go to a minimum 00505 ! we take a guard against too large steps 00506 qs_ot_env(1)%OT_pos(count+1)=MAX(MAXVAL(qs_ot_env(1)%OT_pos(1:3))*0.01_dp, & 00507 MIN(pos,MAXVAL(qs_ot_env(1)%OT_pos(1:3))*4.0_dp)) 00508 ELSE ! just take an extended step 00509 qs_ot_env(1)%OT_pos(count+1)=MAXVAL(qs_ot_env(1)%OT_pos(1:3))*2.0_dp 00510 ENDIF 00511 qs_ot_env(1)%energy_only=.FALSE. 00512 qs_ot_env(1)%line_search_might_be_done=.TRUE. 00513 CASE DEFAULT 00514 CALL stop_program(routineN,moduleN,__LINE__,"NYI") 00515 END SELECT 00516 ds=qs_ot_env(1)%OT_pos(count+1)-qs_ot_env(1)%OT_pos(count) 00517 qs_ot_env(1)%ds_min=qs_ot_env(1)%OT_pos(count+1) 00518 00519 CALL take_step(ds,qs_ot_env,error=error) 00520 00521 END SUBROUTINE do_line_search_3pnt 00522 00523 ! ***************************************************************************** 00524 SUBROUTINE do_line_search_2pnt(qs_ot_env,error) 00525 00526 TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env 00527 TYPE(cp_error_type), INTENT(inout) :: error 00528 00529 CHARACTER(len=*), PARAMETER :: routineN = 'do_line_search_2pnt', 00530 routineP = moduleN//':'//routineN 00531 00532 INTEGER :: count 00533 REAL(KIND=dp) :: a, b, c, ds, pos, val, x0, x1 00534 00535 qs_ot_env(1)%line_search_might_be_done=.FALSE. 00536 qs_ot_env(1)%energy_only=.TRUE. 00537 00538 ! a three point interpolation based on the energy 00539 qs_ot_env(1)%line_search_count=qs_ot_env(1)%line_search_count+1 00540 count=qs_ot_env(1)%line_search_count 00541 qs_ot_env(1)%ot_energy(count)=qs_ot_env(1)%etotal 00542 SELECT CASE(count) 00543 CASE(1) 00544 qs_ot_env(1)%ot_pos(count)=0.0_dp 00545 qs_ot_env(1)%ot_grad(count)=qs_ot_env(1)%gradient 00546 qs_ot_env(1)%ot_pos(count+1)=qs_ot_env(1)%ds_min*1.0_dp 00547 CASE(2) 00548 x0=0.0_dp 00549 c=qs_ot_env(1)%ot_energy(1) 00550 b=qs_ot_env(1)%ot_grad(1) 00551 x1=qs_ot_env(1)%ot_pos(2) 00552 a=(qs_ot_env(1)%ot_energy(2)-b*x1-c)/(x1**2) 00553 IF (a.le.0.0_dp) a=1.0E-15_dp 00554 pos=-b/(2.0_dp*a) 00555 val=a*pos**2+b*pos+c 00556 qs_ot_env(1)%energy_only=.FALSE. 00557 qs_ot_env(1)%line_search_might_be_done=.TRUE. 00558 IF (val.lt.qs_ot_env(1)%ot_energy(1) .AND. val.le.qs_ot_env(1)%ot_energy(2)) THEN 00559 ! we go to a minimum, but ... 00560 ! we take a guard against too large steps 00561 qs_ot_env(1)%OT_pos(count+1)=MAX(MAXVAL(qs_ot_env(1)%OT_pos(1:2))*0.01_dp, & 00562 MIN(pos,MAXVAL(qs_ot_env(1)%OT_pos(1:2))*4.0_dp)) 00563 ELSE ! just take an extended step 00564 qs_ot_env(1)%OT_pos(count+1)=MAXVAL(qs_ot_env(1)%OT_pos(1:2))*2.0_dp 00565 ENDIF 00566 CASE DEFAULT 00567 CALL stop_program(routineN,moduleN,__LINE__,"NYI") 00568 END SELECT 00569 ds=qs_ot_env(1)%OT_pos(count+1)-qs_ot_env(1)%OT_pos(count) 00570 qs_ot_env(1)%ds_min=qs_ot_env(1)%OT_pos(count+1) 00571 00572 CALL take_step(ds,qs_ot_env,error=error) 00573 00574 END SUBROUTINE do_line_search_2pnt 00575 00576 ! ***************************************************************************** 00577 SUBROUTINE do_line_search_none(qs_ot_env,error) 00578 TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env 00579 TYPE(cp_error_type), INTENT(inout) :: error 00580 00581 CALL take_step(qs_ot_env(1)%ds_min,qs_ot_env,error=error) 00582 00583 END SUBROUTINE do_line_search_none 00584 00585 ! 00586 ! creates a new SD direction, using the preconditioner if associated 00587 ! also updates the gradient for line search 00588 ! 00589 00590 ! ***************************************************************************** 00591 SUBROUTINE ot_new_sd_direction(qs_ot_env,error) 00592 TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env 00593 TYPE(cp_error_type), INTENT(inout) :: error 00594 00595 CHARACTER(len=*), PARAMETER :: routineN = 'ot_new_sd_direction', 00596 routineP = moduleN//':'//routineN 00597 00598 INTEGER :: ispin, itmp, k, n, nener, 00599 nscp, nspin 00600 LOGICAL :: do_ener, do_ks, do_scp_dft, 00601 do_scp_nddo 00602 REAL(KIND=dp) :: tmp 00603 TYPE(cp_logger_type), POINTER :: logger 00604 00605 !***SCP 00606 00607 nspin=SIZE(qs_ot_env) 00608 logger=>cp_error_get_logger(error) 00609 do_ks = qs_ot_env ( 1 ) % settings % ks 00610 do_scp_dft = qs_ot_env ( 1 ) % settings % scp_dft 00611 do_scp_nddo = qs_ot_env ( 1 ) % settings % scp_nddo 00612 do_ener = qs_ot_env ( 1 ) % settings % do_ener 00613 00614 IF (ASSOCIATED(qs_ot_env(1)%preconditioner)) THEN 00615 IF (.NOT. qs_ot_env(1)%use_dx) CALL stop_program(routineN,moduleN,__LINE__,"use dx") 00616 qs_ot_env(1)%gnorm=0.0_dp 00617 IF ( do_ks ) THEN 00618 DO ispin=1,nspin 00619 CALL apply_preconditioner(qs_ot_env(ispin)%preconditioner, & 00620 qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_dx,error=error) 00621 CALL cp_dbcsr_trace(qs_ot_env(ispin)%matrix_gx,qs_ot_env(ispin)%matrix_dx,tmp,error=error) 00622 qs_ot_env(1)%gnorm=qs_ot_env(1)%gnorm+tmp 00623 ENDDO 00624 IF (qs_ot_env(1)%gnorm .LT. 0.0_dp) THEN 00625 logger=>cp_error_get_logger(error) 00626 WRITE(cp_logger_get_default_unit_nr(logger),*) "WARNING Preconditioner not positive definite !" 00627 ENDIF 00628 DO ispin=1,nspin 00629 CALL cp_dbcsr_scale(qs_ot_env(ispin)%matrix_dx,-1.0_dp,error=error) 00630 ENDDO 00631 IF (qs_ot_env(1)%settings%do_rotation) THEN 00632 DO ispin=1,nspin 00633 ! right now no preconditioner yet 00634 CALL cp_dbcsr_copy(qs_ot_env(ispin)%rot_mat_dx,qs_ot_env(ispin)%rot_mat_gx,error=error) 00635 CALL cp_dbcsr_trace(qs_ot_env(ispin)%rot_mat_gx,qs_ot_env(ispin)%rot_mat_dx,tmp,error=error) 00636 ! added 0.5, because we have (antisymmetry) only half the number of variables 00637 qs_ot_env(1)%gnorm=qs_ot_env(1)%gnorm+0.5_dp*tmp 00638 ENDDO 00639 DO ispin=1,nspin 00640 CALL cp_dbcsr_scale(qs_ot_env(ispin)%rot_mat_dx,-1.0_dp,error=error) 00641 ENDDO 00642 ENDIF 00643 ENDIF 00644 ! **** SCP 00645 IF (do_scp_dft ) THEN 00646 ! Remember, SCP is formally a spin restricted theory 00647 ! right now no preconditioner yet 00648 qs_ot_env ( 1 ) % dx = qs_ot_env ( 1 ) % gx 00649 tmp = DOT_PRODUCT ( qs_ot_env ( 1 ) % gx, qs_ot_env ( 1 ) % dx ) 00650 qs_ot_env ( 1 ) % dx = -1.0_dp * qs_ot_env ( 1 ) % dx 00651 CALL mp_sum ( tmp, qs_ot_env ( 1 ) % scp_para_env % group ) 00652 qs_ot_env ( 1 )%gnorm=qs_ot_env(1)%gnorm+tmp 00653 ENDIF 00654 IF (do_scp_nddo ) THEN 00655 ! Remember, SCP is formally a spin restricted theory 00656 ! right now no preconditioner yet 00657 CALL cp_dbcsr_copy ( qs_ot_env ( 1 ) % dxmat,qs_ot_env ( 1 ) % gxmat, error=error ) 00658 CALL cp_dbcsr_trace ( qs_ot_env ( 1 ) % gxmat, qs_ot_env ( 1 ) % dxmat, tmp, local_sum=.TRUE.,& 00659 error=error) 00660 CALL cp_dbcsr_scale( qs_ot_env ( 1 ) % dxmat, -1.0_dp , error=error) 00661 qs_ot_env ( 1 )%gnorm=qs_ot_env(1)%gnorm+tmp 00662 ENDIF 00663 ! **** SCP 00664 IF (do_ener) THEN 00665 DO ispin=1,nspin 00666 qs_ot_env(ispin)%ener_dx=qs_ot_env(ispin)%ener_gx 00667 tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_dx,qs_ot_env(ispin)%ener_gx) 00668 qs_ot_env(1)%gnorm=qs_ot_env(1)%gnorm+tmp 00669 qs_ot_env(ispin)%ener_dx=-qs_ot_env(ispin)%ener_dx 00670 ENDDO 00671 ENDIF 00672 ELSE 00673 qs_ot_env(1)%gnorm=0.0_dp 00674 IF ( do_ks ) THEN 00675 DO ispin=1,nspin 00676 CALL cp_dbcsr_trace(qs_ot_env(ispin)%matrix_gx,qs_ot_env(ispin)%matrix_gx,tmp,error=error) 00677 qs_ot_env(1)%gnorm=qs_ot_env(1)%gnorm+tmp 00678 ENDDO 00679 IF (qs_ot_env(1)%settings%do_rotation) THEN 00680 DO ispin=1,nspin 00681 CALL cp_dbcsr_trace(qs_ot_env(ispin)%rot_mat_gx,qs_ot_env(ispin)%rot_mat_gx,tmp,error=error) 00682 ! added 0.5, because we have (antisymmetry) only half the number of variables 00683 qs_ot_env(1)%gnorm=qs_ot_env(1)%gnorm+0.5_dp*tmp 00684 ENDDO 00685 ENDIF 00686 ENDIF 00687 ! **** SCP 00688 IF ( do_scp_dft ) THEN 00689 ! Remember, SCP is formally a spin restricted theory 00690 tmp = DOT_PRODUCT ( qs_ot_env ( 1 ) % gx, qs_ot_env ( 1 ) % gx ) 00691 CALL mp_sum ( tmp, qs_ot_env ( 1 ) % scp_para_env % group ) 00692 qs_ot_env ( 1 )%gnorm=qs_ot_env(1)%gnorm+tmp 00693 ENDIF 00694 IF ( do_scp_nddo ) THEN 00695 ! Remember, SCP is formally a spin restricted theory 00696 CALL cp_dbcsr_trace ( qs_ot_env ( 1 ) % gxmat, qs_ot_env ( 1 ) % gxmat, tmp, local_sum=.TRUE.,& 00697 error=error) 00698 qs_ot_env ( 1 )%gnorm=qs_ot_env(1)%gnorm+tmp 00699 ENDIF 00700 ! **** SCP 00701 IF (do_ener) THEN 00702 DO ispin=1,nspin 00703 tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx,qs_ot_env(ispin)%ener_gx) 00704 qs_ot_env(1)%gnorm=qs_ot_env(1)%gnorm+tmp 00705 ENDDO 00706 ENDIF 00707 ENDIF 00708 00709 k=0 00710 n=0 00711 nscp=0 00712 nener=0 00713 IF ( do_ks ) THEN 00714 CALL cp_dbcsr_get_info(qs_ot_env(1)%matrix_x,nfullrows_total=n) 00715 DO ispin=1,nspin 00716 CALL cp_dbcsr_get_info(qs_ot_env(ispin)%matrix_x,nfullcols_total=itmp) 00717 k=k+itmp 00718 ENDDO 00719 ENDIF 00720 ! **** SCP 00721 IF ( do_scp_nddo ) THEN 00722 nscp = qs_ot_env ( 1 ) % n_el_scp 00723 ELSEIF ( do_scp_dft) THEN 00724 nscp = SIZE ( qs_ot_env ( 1 ) % x ) 00725 CALL mp_sum ( nscp, qs_ot_env ( 1 ) % scp_para_env % group ) 00726 ENDIF 00727 IF ( do_ener ) THEN 00728 DO ispin=1,nspin 00729 nener = nener + SIZE ( qs_ot_env ( ispin ) % ener_x ) 00730 ENDDO 00731 ENDIF 00732 ! **** SCP 00733 ! Handling the case of no free variables to optimize 00734 IF (INT(n,KIND=int_8)*INT(k,KIND=int_8)+nscp+nener /= 0) THEN 00735 qs_ot_env(1)%delta=SQRT(ABS(qs_ot_env(1)%gnorm)/(INT(n,KIND=int_8)*INT(k,KIND=int_8)+nscp+nener)) 00736 qs_ot_env(1)%gradient = - qs_ot_env(1)%gnorm 00737 ELSE 00738 qs_ot_env(1)%delta=0.0_dp 00739 qs_ot_env(1)%gradient=0.0_dp 00740 END IF 00741 END SUBROUTINE ot_new_sd_direction 00742 00743 ! 00744 ! creates a new CG direction. Implements Polak-Ribierre variant 00745 ! using the preconditioner if associated 00746 ! also updates the gradient for line search 00747 ! 00748 ! ***************************************************************************** 00749 SUBROUTINE ot_new_cg_direction(qs_ot_env,error) 00750 TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env 00751 TYPE(cp_error_type), INTENT(inout) :: error 00752 00753 CHARACTER(len=*), PARAMETER :: routineN = 'ot_new_cg_direction', 00754 routineP = moduleN//':'//routineN 00755 00756 INTEGER :: ispin, itmp, k, n, nener, 00757 nscp, nspin 00758 LOGICAL :: do_ener, do_ks, do_scp_dft, 00759 do_scp_nddo 00760 REAL(KIND=dp) :: beta_pr, gnorm_cross, 00761 test_down, tmp 00762 TYPE(cp_logger_type), POINTER :: logger 00763 00764 nspin=SIZE(qs_ot_env) 00765 logger=>cp_error_get_logger(error) 00766 00767 do_ks = qs_ot_env ( 1 ) % settings % ks 00768 do_scp_dft = qs_ot_env ( 1 ) % settings % scp_dft 00769 do_scp_nddo = qs_ot_env ( 1 ) % settings % scp_nddo 00770 do_ener = qs_ot_env ( 1 ) % settings % do_ener 00771 00772 gnorm_cross=0.0_dp 00773 IF ( do_ks ) THEN 00774 DO ispin=1,nspin 00775 CALL cp_dbcsr_trace(qs_ot_env(ispin)%matrix_gx,qs_ot_env(ispin)%matrix_gx_old,tmp,error=error) 00776 gnorm_cross=gnorm_cross+tmp 00777 ENDDO 00778 IF (qs_ot_env(1)%settings%do_rotation) THEN 00779 DO ispin=1,nspin 00780 CALL cp_dbcsr_trace(qs_ot_env(ispin)%rot_mat_gx,qs_ot_env(ispin)%rot_mat_gx_old,tmp,error=error) 00781 ! added 0.5, because we have (antisymmetry) only half the number of variables 00782 gnorm_cross=gnorm_cross+0.5_dp*tmp 00783 ENDDO 00784 ENDIF 00785 END IF 00786 ! ***SCP 00787 IF ( do_scp_dft ) THEN 00788 tmp = DOT_PRODUCT ( qs_ot_env ( 1 ) % gx, qs_ot_env ( 1 ) % gx_old ) 00789 CALL mp_sum ( tmp, qs_ot_env ( 1 ) % scp_para_env % group ) 00790 gnorm_cross = gnorm_cross + tmp 00791 END IF 00792 IF ( do_scp_nddo ) THEN 00793 CALL cp_dbcsr_trace ( qs_ot_env ( 1 ) % gxmat, qs_ot_env ( 1 ) % gx_oldmat, tmp, local_sum=.TRUE.,& 00794 error=error ) 00795 gnorm_cross = gnorm_cross + tmp 00796 END IF 00797 IF ( do_ener ) THEN 00798 DO ispin=1,nspin 00799 tmp = DOT_PRODUCT ( qs_ot_env ( ispin ) % ener_gx, qs_ot_env ( ispin ) % ener_gx_old ) 00800 gnorm_cross = gnorm_cross + tmp 00801 ENDDO 00802 END IF 00803 ! ***SCP 00804 00805 IF (ASSOCIATED(qs_ot_env(1)%preconditioner)) THEN 00806 00807 DO ispin=1,nspin 00808 CALL apply_preconditioner(qs_ot_env(ispin)%preconditioner, & 00809 qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx_old,error=error) 00810 ENDDO 00811 qs_ot_env(1)%gnorm=0.0_dp 00812 IF ( do_ks ) THEN 00813 DO ispin=1,nspin 00814 CALL cp_dbcsr_trace(qs_ot_env(ispin)%matrix_gx,qs_ot_env(ispin)%matrix_gx_old,tmp,error=error) 00815 qs_ot_env(1)%gnorm=qs_ot_env(1)%gnorm+tmp 00816 ENDDO 00817 IF (qs_ot_env(1)%gnorm .LT. 0.0_dp) THEN 00818 WRITE(cp_logger_get_default_unit_nr(logger),*) "WARNING Preconditioner not positive definite !" 00819 ENDIF 00820 DO ispin=1,nspin 00821 CALL cp_dbcsr_copy(qs_ot_env(ispin)%matrix_gx,qs_ot_env(ispin)%matrix_gx_old,error=error) 00822 ENDDO 00823 IF (qs_ot_env(1)%settings%do_rotation) THEN 00824 DO ispin=1,nspin 00825 ! right now no preconditioner yet 00826 CALL cp_dbcsr_copy(qs_ot_env(ispin)%rot_mat_gx_old,qs_ot_env(ispin)%rot_mat_gx,error=error) 00827 CALL cp_dbcsr_trace(qs_ot_env(ispin)%rot_mat_gx,qs_ot_env(ispin)%rot_mat_gx_old,tmp,error=error) 00828 ! added 0.5, because we have (antisymmetry) only half the number of variables 00829 qs_ot_env(1)%gnorm=qs_ot_env(1)%gnorm+0.5_dp*tmp 00830 ENDDO 00831 DO ispin=1,nspin 00832 CALL cp_dbcsr_copy(qs_ot_env(ispin)%rot_mat_gx,qs_ot_env(ispin)%rot_mat_gx_old,error=error) 00833 ENDDO 00834 ENDIF 00835 END IF 00836 ! ***SCP 00837 IF ( do_scp_dft ) THEN 00838 qs_ot_env ( 1 ) % gx_old = qs_ot_env ( 1 ) % gx 00839 tmp = DOT_PRODUCT ( qs_ot_env ( 1 ) % gx, qs_ot_env ( 1 ) % gx_old ) 00840 CALL mp_sum ( tmp, qs_ot_env ( 1 ) % scp_para_env % group ) 00841 qs_ot_env ( 1 ) % gnorm = qs_ot_env ( 1 ) % gnorm + tmp 00842 qs_ot_env ( 1 ) % gx = qs_ot_env ( 1 ) % gx_old 00843 END IF 00844 IF ( do_scp_nddo ) THEN 00845 CALL cp_dbcsr_copy ( qs_ot_env ( 1 ) % gx_oldmat, qs_ot_env ( 1 ) % gxmat, error=error ) 00846 CALL cp_dbcsr_trace ( qs_ot_env ( 1 ) % gxmat, qs_ot_env ( 1 ) % gx_oldmat, tmp, local_sum=.TRUE.,& 00847 error=error ) 00848 qs_ot_env ( 1 ) % gnorm = qs_ot_env ( 1 ) % gnorm + tmp 00849 CALL cp_dbcsr_copy ( qs_ot_env ( 1 ) % gxmat, qs_ot_env ( 1 ) % gx_oldmat,error=error ) 00850 END IF 00851 IF ( do_ener ) THEN 00852 DO ispin=1,nspin 00853 qs_ot_env ( ispin ) % ener_gx_old = qs_ot_env ( ispin ) % ener_gx 00854 tmp = DOT_PRODUCT ( qs_ot_env ( ispin ) % ener_gx, qs_ot_env ( ispin ) % ener_gx_old ) 00855 qs_ot_env ( 1 ) % gnorm = qs_ot_env ( 1 ) % gnorm + tmp 00856 qs_ot_env ( ispin ) % ener_gx = qs_ot_env ( ispin ) % ener_gx_old 00857 ENDDO 00858 END IF 00859 ELSE 00860 IF ( do_ks ) THEN 00861 qs_ot_env(1)%gnorm=0.0_dp 00862 DO ispin=1,nspin 00863 CALL cp_dbcsr_trace(qs_ot_env(ispin)%matrix_gx,qs_ot_env(ispin)%matrix_gx,tmp,error=error) 00864 qs_ot_env(1)%gnorm=qs_ot_env(1)%gnorm+tmp 00865 CALL cp_dbcsr_copy(qs_ot_env(ispin)%matrix_gx_old,qs_ot_env(ispin)%matrix_gx,error=error) 00866 ENDDO 00867 IF (qs_ot_env(1)%settings%do_rotation) THEN 00868 DO ispin=1,nspin 00869 CALL cp_dbcsr_trace(qs_ot_env(ispin)%rot_mat_gx,qs_ot_env(ispin)%rot_mat_gx,tmp,error=error) 00870 ! added 0.5, because we have (antisymmetry) only half the number of variables 00871 qs_ot_env(1)%gnorm=qs_ot_env(1)%gnorm+0.5_dp*tmp 00872 CALL cp_dbcsr_copy(qs_ot_env(ispin)%rot_mat_gx_old,qs_ot_env(ispin)%rot_mat_gx,error=error) 00873 ENDDO 00874 ENDIF 00875 ENDIF 00876 ! ***SCP 00877 IF ( do_scp_dft ) THEN 00878 tmp = DOT_PRODUCT ( qs_ot_env ( 1 ) % gx, qs_ot_env ( 1 ) % gx ) 00879 CALL mp_sum ( tmp, qs_ot_env ( 1 ) % scp_para_env % group ) 00880 qs_ot_env ( 1 ) % gnorm = qs_ot_env ( 1 ) % gnorm + tmp 00881 qs_ot_env ( 1 ) % gx_old = qs_ot_env ( 1 ) % gx 00882 END IF 00883 IF ( do_scp_nddo ) THEN 00884 CALL cp_dbcsr_trace ( qs_ot_env ( 1 ) % gxmat, qs_ot_env ( 1 ) % gxmat, tmp, local_sum=.TRUE., & 00885 error=error ) 00886 qs_ot_env ( 1 ) % gnorm = qs_ot_env ( 1 ) % gnorm + tmp 00887 CALL cp_dbcsr_copy ( qs_ot_env ( 1 ) % gx_oldmat, qs_ot_env ( 1 ) % gxmat, error=error ) 00888 END IF 00889 ! ***SCP 00890 IF (do_ener) THEN 00891 DO ispin=1,nspin 00892 tmp = DOT_PRODUCT ( qs_ot_env ( ispin ) % ener_gx, qs_ot_env ( ispin ) % ener_gx ) 00893 qs_ot_env ( 1 ) % gnorm = qs_ot_env ( 1 ) % gnorm + tmp 00894 qs_ot_env ( ispin ) % ener_gx_old = qs_ot_env ( ispin ) % ener_gx 00895 ENDDO 00896 ENDIF 00897 ENDIF 00898 00899 k=0 00900 n=0 00901 nscp=0 00902 nener=0 00903 IF ( do_ks ) THEN 00904 CALL cp_dbcsr_get_info(qs_ot_env(1)%matrix_x,nfullrows_total=n) 00905 DO ispin=1,nspin 00906 CALL cp_dbcsr_get_info(qs_ot_env(ispin)%matrix_x,nfullcols_total=itmp) 00907 k=k+itmp 00908 ENDDO 00909 END IF 00910 !***SCP 00911 IF ( do_scp_dft ) THEN 00912 nscp = SIZE ( qs_ot_env ( 1 ) % x ) 00913 CALL mp_sum ( nscp, qs_ot_env ( 1 ) % scp_para_env % group ) 00914 ELSEIF ( do_scp_nddo ) THEN 00915 nscp = qs_ot_env ( 1 ) % n_el_scp 00916 ENDIF 00917 !***SCP 00918 IF (do_ener) THEN 00919 DO ispin=1,nspin 00920 nener=nener+SIZE(qs_ot_env ( ispin ) % ener_x) 00921 ENDDO 00922 ENDIF 00923 ! Handling the case of no free variables to optimize 00924 IF (INT(n,KIND=int_8)*INT(k,KIND=int_8)+nscp+nener /= 0) THEN 00925 qs_ot_env(1)%delta=SQRT(ABS(qs_ot_env(1)%gnorm)/(INT(n,KIND=int_8)*INT(k,KIND=int_8)+nscp+nener)) 00926 beta_pr=(qs_ot_env(1)%gnorm-gnorm_cross)/qs_ot_env(1)%gnorm_old 00927 ELSE 00928 qs_ot_env(1)%delta=0.0_dp 00929 beta_pr = 0.0_dp 00930 END IF 00931 beta_pr=MAX(beta_pr,0.0_dp) ! reset to SD 00932 00933 test_down=0.0_dp 00934 IF ( do_ks ) THEN 00935 DO ispin=1,nspin 00936 CALL cp_dbcsr_add(qs_ot_env(ispin)%matrix_dx,qs_ot_env(ispin)%matrix_gx,& 00937 alpha_scalar=beta_pr,beta_scalar=-1.0_dp,error=error) 00938 CALL cp_dbcsr_trace(qs_ot_env(ispin)%matrix_gx,qs_ot_env(ispin)%matrix_dx,tmp,error=error) 00939 test_down=test_down+tmp 00940 IF (qs_ot_env(1)%settings%do_rotation) THEN 00941 CALL cp_dbcsr_add(qs_ot_env(ispin)%rot_mat_dx,qs_ot_env(ispin)%rot_mat_gx,& 00942 alpha_scalar=beta_pr,beta_scalar=-1.0_dp,error=error) 00943 CALL cp_dbcsr_trace(qs_ot_env(ispin)%rot_mat_gx,qs_ot_env(ispin)%rot_mat_dx,tmp,error=error) 00944 test_down=test_down+0.5_dp*tmp 00945 ENDIF 00946 ENDDO 00947 END IF 00948 ! ***SCP 00949 IF ( do_scp_dft ) THEN 00950 qs_ot_env ( 1 ) % dx = beta_pr * qs_ot_env ( 1 ) % dx - qs_ot_env ( 1 ) % gx 00951 tmp = DOT_PRODUCT ( qs_ot_env ( 1 ) % gx, qs_ot_env ( 1 ) % dx ) 00952 CALL mp_sum ( tmp, qs_ot_env ( 1 ) % scp_para_env % group ) 00953 test_down = test_down + tmp 00954 END IF 00955 IF ( do_scp_nddo ) THEN 00956 CALL cp_dbcsr_add ( qs_ot_env ( 1 ) % dxmat, qs_ot_env ( 1 ) % gxmat, alpha_scalar=beta_pr,& 00957 beta_scalar=-1.0_dp, error=error ) 00958 CALL cp_dbcsr_trace ( qs_ot_env ( 1 ) % gxmat, qs_ot_env ( 1 ) % dxmat, tmp, local_sum=.TRUE., & 00959 error=error ) 00960 test_down = test_down + tmp 00961 END IF 00962 ! ***SCP 00963 IF (do_ener) THEN 00964 DO ispin=1,nspin 00965 qs_ot_env ( ispin ) % ener_dx = beta_pr * qs_ot_env ( ispin ) % ener_dx - qs_ot_env ( ispin ) % ener_gx 00966 tmp = DOT_PRODUCT ( qs_ot_env ( ispin ) % ener_gx, qs_ot_env ( ispin ) % ener_dx ) 00967 test_down = test_down + tmp 00968 ENDDO 00969 ENDIF 00970 00971 IF (test_down.ge.0.0_dp) THEN ! reset to SD 00972 beta_pr=0.0_dp 00973 IF ( do_ks ) THEN 00974 DO ispin=1,nspin 00975 CALL cp_dbcsr_add(qs_ot_env(ispin)%matrix_dx,qs_ot_env(ispin)%matrix_gx,& 00976 alpha_scalar=beta_pr,beta_scalar=-1.0_dp,error=error) 00977 IF (qs_ot_env(1)%settings%do_rotation) THEN 00978 CALL cp_dbcsr_add(qs_ot_env(ispin)%rot_mat_dx, & 00979 qs_ot_env(ispin)%rot_mat_gx,alpha_scalar=beta_pr,beta_scalar=-1.0_dp,error=error) 00980 ENDIF 00981 ENDDO 00982 END IF 00983 ! ***SCP 00984 IF ( do_scp_dft ) THEN 00985 qs_ot_env ( 1 ) % dx = beta_pr * qs_ot_env ( 1 ) % dx - qs_ot_env ( 1 ) % gx 00986 END IF 00987 IF ( do_scp_nddo ) THEN 00988 CALL cp_dbcsr_add ( qs_ot_env ( 1 ) % dxmat, qs_ot_env ( 1 ) % gxmat, & 00989 alpha_scalar=beta_pr, beta_scalar=-1.0_dp, error=error ) 00990 END IF 00991 ! ***SCP 00992 IF ( do_ener ) THEN 00993 DO ispin=1,nspin 00994 qs_ot_env ( ispin ) % ener_dx = beta_pr * qs_ot_env ( ispin ) % ener_dx - qs_ot_env ( ispin ) % ener_gx 00995 ENDDO 00996 ENDIF 00997 ENDIF 00998 ! since we change the direction we have to adjust the gradient 00999 qs_ot_env(1)%gradient = beta_pr*qs_ot_env(1)%gradient - qs_ot_env(1)%gnorm 01000 qs_ot_env(1)%gnorm_old=qs_ot_env(1)%gnorm 01001 END SUBROUTINE ot_new_cg_direction 01002 01003 ! ***************************************************************************** 01004 SUBROUTINE ot_diis_step(qs_ot_env,error) 01005 TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env 01006 TYPE(cp_error_type), INTENT(inout) :: error 01007 01008 CHARACTER(len=*), PARAMETER :: routineN = 'ot_diis_step', 01009 routineP = moduleN//':'//routineN 01010 01011 INTEGER :: diis_bound, diis_m, i, info, 01012 ispin, itmp, j, k, n, nener, 01013 nscp, nspin 01014 LOGICAL :: do_ener, do_ks, do_scp_dft, 01015 do_scp_nddo 01016 REAL(KIND=dp) :: overlap, tmp, tr_xnew_gx, 01017 tr_xold_gx 01018 TYPE(cp_logger_type), POINTER :: logger 01019 01020 do_ks = qs_ot_env ( 1 ) % settings % ks 01021 do_scp_dft = qs_ot_env ( 1 ) % settings % scp_dft 01022 do_scp_nddo = qs_ot_env ( 1 ) % settings % scp_nddo 01023 do_ener = qs_ot_env ( 1 ) % settings % do_ener 01024 nspin=SIZE(qs_ot_env) 01025 01026 diis_m=qs_ot_env(1)%settings%diis_m 01027 01028 IF (qs_ot_env(1)%diis_iter.lt.diis_m) THEN 01029 diis_bound=qs_ot_env(1)%diis_iter+1 01030 ELSE 01031 diis_bound=diis_m 01032 ENDIF 01033 01034 j = MOD(qs_ot_env(1)%diis_iter,diis_m)+1 ! index in the circular array 01035 01036 ! copy the position and the error vector in the diis buffers 01037 01038 IF ( do_ks ) THEN 01039 DO ispin=1,nspin 01040 CALL cp_dbcsr_copy(qs_ot_env(ispin)%matrix_h_x(j)%matrix,qs_ot_env(ispin)%matrix_x,error=error) 01041 IF (qs_ot_env(ispin)%settings%do_rotation) THEN 01042 CALL cp_dbcsr_copy(qs_ot_env(ispin)%rot_mat_h_x(j)%matrix,qs_ot_env(ispin)%rot_mat_x,error=error) 01043 ENDIF 01044 ENDDO 01045 END IF 01046 IF ( do_scp_dft )THEN 01047 qs_ot_env ( 1 ) % h_x ( j, : ) = qs_ot_env ( 1 ) % x ( : ) 01048 END IF 01049 IF ( do_scp_nddo ) THEN 01050 CALL cp_dbcsr_copy ( qs_ot_env ( 1 ) % hx_mat ( j )%matrix, qs_ot_env ( 1 ) % xmat, error=error) 01051 ENDIF 01052 IF ( do_ener )THEN 01053 DO ispin=1,nspin 01054 qs_ot_env ( ispin ) % ener_h_x ( j, : ) = qs_ot_env ( ispin ) % ener_x ( : ) 01055 ENDDO 01056 END IF 01057 ! *** SCP 01058 IF (ASSOCIATED(qs_ot_env(1)%preconditioner)) THEN 01059 qs_ot_env(1)%gnorm=0.0_dp 01060 IF ( do_ks ) THEN 01061 DO ispin=1,nspin 01062 CALL apply_preconditioner(qs_ot_env(ispin)%preconditioner, & 01063 qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_h_e(j)%matrix,error=error) 01064 CALL cp_dbcsr_trace(qs_ot_env(ispin)%matrix_gx,qs_ot_env(ispin)%matrix_h_e(j)%matrix, & 01065 tmp,error=error) 01066 qs_ot_env(1)%gnorm=qs_ot_env(1)%gnorm+tmp 01067 ENDDO 01068 IF (qs_ot_env(1)%gnorm .LT. 0.0_dp) THEN 01069 WRITE(cp_logger_get_default_unit_nr(logger),*) "WARNING Preconditioner not positive definite !" 01070 ENDIF 01071 DO ispin=1,nspin 01072 CALL cp_dbcsr_scale(qs_ot_env(ispin)%matrix_h_e(j)%matrix,-qs_ot_env(1)%ds_min,error=error) 01073 ENDDO 01074 IF (qs_ot_env(1)%settings%do_rotation) THEN 01075 DO ispin=1,nspin 01076 CALL cp_dbcsr_copy(qs_ot_env(ispin)%rot_mat_h_e(j)%matrix,qs_ot_env(ispin)%rot_mat_gx,error=error) 01077 CALL cp_dbcsr_trace(qs_ot_env(ispin)%rot_mat_gx,qs_ot_env(ispin)%rot_mat_h_e(j)%matrix, & 01078 tmp,error=error) 01079 qs_ot_env(1)%gnorm=qs_ot_env(1)%gnorm+0.5_dp*tmp 01080 ENDDO 01081 DO ispin=1,nspin 01082 CALL cp_dbcsr_scale(qs_ot_env(ispin)%rot_mat_h_e(j)%matrix,-qs_ot_env(1)%ds_min,error=error) 01083 ENDDO 01084 ENDIF 01085 END IF 01086 ! ***SCP 01087 IF (do_scp_dft) THEN 01088 qs_ot_env ( 1 ) % h_e ( j, : ) = qs_ot_env ( 1 ) % gx ( : ) 01089 tmp = DOT_PRODUCT ( qs_ot_env ( 1 ) % h_e ( j, : ), qs_ot_env ( 1 ) % gx ( : ) ) 01090 CALL mp_sum ( tmp, qs_ot_env ( 1 ) % scp_para_env % group ) 01091 qs_ot_env ( 1 ) % gnorm = qs_ot_env ( 1 ) % gnorm + tmp 01092 qs_ot_env ( 1 ) % h_e ( j, : ) = -qs_ot_env ( 1 ) % ds_min * qs_ot_env ( 1 ) % h_e ( j, : ) 01093 END IF 01094 IF ( do_scp_nddo ) THEN 01095 CALL cp_dbcsr_copy ( qs_ot_env ( 1 ) % he_mat ( j ) % matrix ,qs_ot_env ( 1 ) % gxmat,error=error) 01096 CALL cp_dbcsr_trace ( qs_ot_env ( 1 ) % gxmat, qs_ot_env ( 1 ) % he_mat ( j ) % matrix, tmp, & 01097 local_sum=.TRUE. ,error=error ) 01098 qs_ot_env ( 1 ) % gnorm = qs_ot_env ( 1 ) % gnorm + tmp 01099 CALL cp_dbcsr_scale ( qs_ot_env ( 1 ) % he_mat ( j ) % matrix, -qs_ot_env ( 1 ) % ds_min, error=error ) 01100 ENDIF 01101 ! ***SCP 01102 IF (do_ener) THEN 01103 DO ispin=1,nspin 01104 qs_ot_env ( ispin ) % ener_h_e ( j, : ) = qs_ot_env ( ispin ) % ener_gx ( : ) 01105 tmp = DOT_PRODUCT ( qs_ot_env ( ispin ) % ener_h_e ( j, : ), qs_ot_env ( ispin ) % ener_gx ( : ) ) 01106 qs_ot_env ( 1 ) % gnorm = qs_ot_env ( 1 ) % gnorm + tmp 01107 qs_ot_env ( ispin ) % ener_h_e ( j, : ) = -qs_ot_env ( 1 ) % ds_min * qs_ot_env ( ispin ) % ener_h_e ( j, : ) 01108 ENDDO 01109 ENDIF 01110 ELSE 01111 qs_ot_env(1)%gnorm=0.0_dp 01112 IF ( do_ks ) THEN 01113 DO ispin=1,nspin 01114 CALL cp_dbcsr_trace(qs_ot_env(ispin)%matrix_gx,qs_ot_env(ispin)%matrix_gx,tmp,error=error) 01115 qs_ot_env(1)%gnorm=qs_ot_env(1)%gnorm+tmp 01116 CALL cp_dbcsr_add(qs_ot_env(ispin)%matrix_h_e(j)%matrix, & 01117 qs_ot_env(ispin)%matrix_gx,alpha_scalar=0.0_dp,beta_scalar=-qs_ot_env(1)%ds_min,error=error) 01118 ENDDO 01119 IF (qs_ot_env(1)%settings%do_rotation) THEN 01120 DO ispin=1,nspin 01121 CALL cp_dbcsr_trace(qs_ot_env(ispin)%rot_mat_gx,qs_ot_env(ispin)%rot_mat_gx,tmp,error=error) 01122 qs_ot_env(1)%gnorm=qs_ot_env(1)%gnorm+0.5_dp*tmp 01123 CALL cp_dbcsr_add(qs_ot_env(ispin)%rot_mat_h_e(j)%matrix, & 01124 qs_ot_env(ispin)%rot_mat_gx,alpha_scalar=0.0_dp,beta_scalar=-qs_ot_env(1)%ds_min,error=error) 01125 ENDDO 01126 ENDIF 01127 END IF 01128 ! ***SCP 01129 IF (do_scp_dft) THEN 01130 tmp = DOT_PRODUCT ( qs_ot_env ( 1 ) % gx ( : ), qs_ot_env ( 1 ) % gx ( : ) ) 01131 CALL mp_sum ( tmp, qs_ot_env ( 1 ) % scp_para_env % group ) 01132 qs_ot_env ( 1 ) % gnorm = qs_ot_env ( 1 ) % gnorm + tmp 01133 qs_ot_env ( 1 ) % h_e ( j, : ) = -qs_ot_env ( 1 ) % ds_min * qs_ot_env ( 1 ) % gx ( : ) 01134 END IF 01135 IF ( do_scp_nddo ) THEN 01136 CALL cp_dbcsr_trace ( qs_ot_env ( 1 ) % gxmat, qs_ot_env ( 1 ) % gxmat, tmp, local_sum=.TRUE., & 01137 error=error ) 01138 qs_ot_env ( 1 ) % gnorm = qs_ot_env ( 1 ) % gnorm + tmp 01139 CALL cp_dbcsr_copy ( qs_ot_env ( 1 ) % he_mat ( j )%matrix,qs_ot_env ( 1 ) % gxmat,error=error) 01140 CALL cp_dbcsr_add ( qs_ot_env ( 1 ) % he_mat ( j ) % matrix, & 01141 qs_ot_env ( 1 ) % gxmat,alpha_scalar=0.0_dp,beta_scalar=-qs_ot_env ( 1 ) % ds_min, error=error ) 01142 END IF 01143 ! ***SCP 01144 IF (do_ener) THEN 01145 DO ispin=1,nspin 01146 tmp = DOT_PRODUCT ( qs_ot_env ( ispin ) % ener_gx ( : ), qs_ot_env ( ispin ) % ener_gx ( : ) ) 01147 qs_ot_env ( 1 ) % gnorm = qs_ot_env ( 1 ) % gnorm + tmp 01148 qs_ot_env ( ispin ) % ener_h_e ( j, : ) = -qs_ot_env ( 1 ) % ds_min * qs_ot_env ( ispin ) % ener_gx ( : ) 01149 ENDDO 01150 END IF 01151 ENDIF 01152 k = 0 01153 nscp = 0 01154 n = 0 01155 nener= 0 01156 IF ( do_ks ) THEN 01157 CALL cp_dbcsr_get_info(qs_ot_env(1)%matrix_x,nfullrows_total=n) 01158 DO ispin=1,nspin 01159 CALL cp_dbcsr_get_info(qs_ot_env(ispin)%matrix_x,nfullcols_total=itmp) 01160 k=k+itmp 01161 ENDDO 01162 END IF 01163 !***SCP 01164 IF ( do_scp_dft ) THEN 01165 nscp = SIZE ( qs_ot_env ( 1 ) % x ) 01166 CALL mp_sum ( nscp, qs_ot_env ( 1 ) % scp_para_env % group ) 01167 ELSEIF ( do_scp_nddo ) THEN 01168 nscp = qs_ot_env ( 1 ) % n_el_scp 01169 ENDIF 01170 !***SCP 01171 IF ( do_ener ) THEN 01172 DO ispin=1,nspin 01173 nener = nener + SIZE( qs_ot_env ( ispin ) % ener_x ) 01174 ENDDO 01175 ENDIF 01176 ! Handling the case of no free variables to optimize 01177 IF (INT(n,KIND=int_8)*INT(k,KIND=int_8)+nscp+nener /= 0) THEN 01178 qs_ot_env(1)%delta=SQRT(ABS(qs_ot_env(1)%gnorm)/(INT(n,KIND=int_8)*INT(k,KIND=int_8)+nscp+nener)) 01179 qs_ot_env(1)%gradient = - qs_ot_env(1)%gnorm 01180 ELSE 01181 qs_ot_env(1)%delta = 0.0_dp 01182 qs_ot_env(1)%gradient = 0.0_dp 01183 END IF 01184 01185 ! make the diis matrix and solve it 01186 DO i=1,diis_bound 01187 ! I think there are two possible options, with and without preconditioner 01188 ! as a metric 01189 ! the second option seems most logical to me, and it seems marginally faster 01190 ! in some of the tests 01191 IF (.FALSE.) THEN 01192 qs_ot_env(1)%ls_diis(i,j)=0.0_dp 01193 IF ( do_ks ) THEN 01194 DO ispin=1,nspin 01195 CALL cp_dbcsr_trace(qs_ot_env(ispin)%matrix_h_e(j)%matrix, & 01196 qs_ot_env(ispin)%matrix_h_e(i)%matrix, & 01197 tmp,error=error) 01198 qs_ot_env(1)%ls_diis(i,j)=qs_ot_env(1)%ls_diis(i,j)+tmp 01199 IF (qs_ot_env(ispin)%settings%do_rotation) THEN 01200 CALL cp_dbcsr_trace(qs_ot_env(ispin)%rot_mat_h_e(j)%matrix, & 01201 qs_ot_env(ispin)%rot_mat_h_e(i)%matrix, & 01202 tmp,error=error) 01203 qs_ot_env(1)%ls_diis(i,j)=qs_ot_env(1)%ls_diis(i,j)+0.5_dp*tmp 01204 ENDIF 01205 ENDDO 01206 END IF 01207 ! ***SCP 01208 IF (do_scp_dft) THEN 01209 tmp = DOT_PRODUCT ( qs_ot_env ( 1 ) % h_e ( j, : ), qs_ot_env ( 1 ) % h_e ( i, : ) ) 01210 CALL mp_sum ( tmp, qs_ot_env ( 1 ) % scp_para_env % group ) 01211 qs_ot_env ( 1 ) % ls_diis ( i, j ) = qs_ot_env ( 1 ) % ls_diis ( i, j ) + tmp 01212 END IF 01213 IF ( do_scp_nddo ) THEN 01214 CALL cp_dbcsr_trace ( qs_ot_env ( 1 ) % he_mat ( j ) % matrix, qs_ot_env ( 1 ) % he_mat ( i ) % matrix, & 01215 tmp, local_sum=.TRUE.,error=error ) 01216 qs_ot_env ( 1 ) % ls_diis ( i, j ) = qs_ot_env ( 1 ) % ls_diis ( i, j ) + tmp 01217 END IF 01218 IF (do_ener) THEN 01219 DO ispin=1,nspin 01220 tmp = DOT_PRODUCT ( qs_ot_env ( ispin ) % ener_h_e ( j, : ), qs_ot_env ( ispin ) % ener_h_e ( i, : ) ) 01221 qs_ot_env ( 1 ) % ls_diis ( i, j ) = qs_ot_env ( 1 ) % ls_diis ( i, j ) + tmp 01222 ENDDO 01223 END IF 01224 ! ***SCP 01225 ELSE 01226 qs_ot_env(1)%ls_diis(i,j)=0.0_dp 01227 IF ( do_ks ) THEN 01228 DO ispin=1,nspin 01229 CALL cp_dbcsr_trace(qs_ot_env(ispin)%matrix_gx, & 01230 qs_ot_env(ispin)%matrix_h_e(i)%matrix, & 01231 tmp,error=error) 01232 qs_ot_env(1)%ls_diis(i,j)=qs_ot_env(1)%ls_diis(i,j)-qs_ot_env(1)%ds_min * tmp 01233 IF (qs_ot_env(ispin)%settings%do_rotation) THEN 01234 CALL cp_dbcsr_trace(qs_ot_env(ispin)%rot_mat_gx, & 01235 qs_ot_env(ispin)%rot_mat_h_e(i)%matrix, & 01236 tmp,error=error) 01237 qs_ot_env(1)%ls_diis(i,j)=qs_ot_env(1)%ls_diis(i,j)-qs_ot_env(1)%ds_min * 0.5_dp * tmp 01238 ENDIF 01239 ENDDO 01240 END IF 01241 ! ***SCP 01242 IF (do_scp_dft) THEN 01243 tmp = DOT_PRODUCT ( qs_ot_env ( 1 ) % gx ( : ), qs_ot_env ( 1 ) % h_e ( i, : ) ) 01244 CALL mp_sum ( tmp, qs_ot_env ( 1 ) % scp_para_env % group ) 01245 qs_ot_env ( 1 ) % ls_diis ( i, j ) = qs_ot_env ( 1 ) % ls_diis ( i, j ) - qs_ot_env ( 1 ) % ds_min * tmp 01246 END IF 01247 IF ( do_scp_nddo ) THEN 01248 CALL cp_dbcsr_trace ( qs_ot_env ( 1 ) % gxmat, qs_ot_env ( 1 ) % he_mat ( i ) % matrix, tmp, & 01249 local_sum=.TRUE., error=error ) 01250 qs_ot_env ( 1 ) % ls_diis ( i, j ) = qs_ot_env ( 1 ) % ls_diis ( i, j ) - qs_ot_env ( 1 ) % ds_min * tmp 01251 END IF 01252 ! ***SCP 01253 IF (do_ener) THEN 01254 DO ispin=1,nspin 01255 tmp = DOT_PRODUCT ( qs_ot_env ( ispin ) % ener_gx ( : ), qs_ot_env ( ispin ) % ener_h_e ( i, : ) ) 01256 qs_ot_env ( 1 ) % ls_diis ( i, j ) = qs_ot_env ( 1 ) % ls_diis ( i, j ) - qs_ot_env ( 1 )% ds_min * tmp 01257 ENDDO 01258 END IF 01259 ENDIF 01260 qs_ot_env(1)%ls_diis(j,i)=qs_ot_env(1)%ls_diis(i,j) 01261 qs_ot_env(1)%ls_diis(i,diis_bound+1)=1.0_dp 01262 qs_ot_env(1)%ls_diis(diis_bound+1,i)=1.0_dp 01263 qs_ot_env(1)%c_diis(i)=0.0_dp 01264 ENDDO 01265 qs_ot_env(1)%ls_diis(diis_bound+1,diis_bound+1)=0.0_dp 01266 qs_ot_env(1)%c_diis(diis_bound+1)=1.0_dp 01267 ! put in buffer, dgesv destroys 01268 qs_ot_env(1)%lss_diis=qs_ot_env(1)%ls_diis 01269 01270 CALL DGESV(diis_bound+1, 1, qs_ot_env(1)%lss_diis,diis_m+1,qs_ot_env(1)%ipivot,& 01271 qs_ot_env(1)%c_diis, diis_m+1, info) 01272 IF (info.ne.0) CALL stop_program(routineN,moduleN,__LINE__,"Singular DIIS matrix") 01273 01274 IF ( do_ks ) THEN 01275 DO ispin=1,nspin 01276 ! OK, add the vectors now 01277 CALL cp_dbcsr_set(qs_ot_env(ispin)%matrix_x,0.0_dp,error=error) 01278 DO i=1, diis_bound 01279 CALL cp_dbcsr_add(qs_ot_env(ispin)%matrix_x, & 01280 qs_ot_env(ispin)%matrix_h_e(i)%matrix,& 01281 alpha_scalar=1.0_dp,beta_scalar=qs_ot_env(1)%c_diis(i),error=error) 01282 ENDDO 01283 DO i=1, diis_bound 01284 CALL cp_dbcsr_add(qs_ot_env(ispin)%matrix_x, & 01285 qs_ot_env(ispin)%matrix_h_x(i)%matrix,& 01286 alpha_scalar=1.0_dp,beta_scalar=qs_ot_env(1)%c_diis(i),error=error) 01287 ENDDO 01288 IF (qs_ot_env(ispin)%settings%do_rotation) THEN 01289 CALL cp_dbcsr_set(qs_ot_env(ispin)%rot_mat_x,0.0_dp,error=error) 01290 DO i=1, diis_bound 01291 CALL cp_dbcsr_add(qs_ot_env(ispin)%rot_mat_x, & 01292 qs_ot_env(ispin)%rot_mat_h_e(i)%matrix,& 01293 alpha_scalar=1.0_dp,beta_scalar=qs_ot_env(1)%c_diis(i),error=error) 01294 ENDDO 01295 DO i=1, diis_bound 01296 CALL cp_dbcsr_add(qs_ot_env(ispin)%rot_mat_x, & 01297 qs_ot_env(ispin)%rot_mat_h_x(i)%matrix,& 01298 alpha_scalar=1.0_dp,beta_scalar=qs_ot_env(1)%c_diis(i),error=error) 01299 ENDDO 01300 ENDIF 01301 ENDDO 01302 END IF 01303 ! ***SCP 01304 IF (do_scp_dft) THEN 01305 qs_ot_env ( 1 ) % x ( : ) = 0.0_dp 01306 DO i = 1, diis_bound 01307 qs_ot_env ( 1 ) % x ( : ) = qs_ot_env ( 1 ) % x ( : ) & 01308 + qs_ot_env ( 1 ) % c_diis ( i ) * qs_ot_env ( 1 ) % h_e ( i, : ) 01309 END DO 01310 DO i = 1, diis_bound 01311 qs_ot_env ( 1 ) % x ( : ) = qs_ot_env ( 1 ) % x ( : ) & 01312 + qs_ot_env ( 1 ) % c_diis ( i ) * qs_ot_env ( 1 ) % h_x ( i, : ) 01313 END DO 01314 END IF 01315 IF (do_scp_nddo) THEN 01316 CALL cp_dbcsr_set( qs_ot_env ( 1 ) % xmat, 0.0_dp, error=error) 01317 DO i = 1, diis_bound 01318 CALL cp_dbcsr_add ( qs_ot_env ( 1 ) % xmat, qs_ot_env ( 1 ) % he_mat ( i ) % matrix, & 01319 alpha_scalar=1.0_dp, beta_scalar=qs_ot_env ( 1 ) % c_diis ( i ), error=error ) 01320 END DO 01321 DO i = 1, diis_bound 01322 CALL cp_dbcsr_add ( qs_ot_env ( 1 ) % xmat, qs_ot_env ( 1 ) % hx_mat ( i ) % matrix, & 01323 alpha_scalar=1.0_dp,beta_scalar=qs_ot_env ( 1 ) % c_diis ( i ), error=error ) 01324 END DO 01325 END IF 01326 ! ***SCP 01327 IF (do_ener) THEN 01328 DO ispin=1,nspin 01329 qs_ot_env ( ispin ) % ener_x ( : ) = 0.0_dp 01330 DO i = 1, diis_bound 01331 qs_ot_env ( ispin ) % ener_x ( : ) = qs_ot_env ( ispin ) % ener_x ( : ) & 01332 + qs_ot_env ( 1 ) % c_diis ( i ) * qs_ot_env ( ispin ) % ener_h_e ( i, : ) 01333 END DO 01334 DO i = 1, diis_bound 01335 qs_ot_env ( ispin ) % ener_x ( : ) = qs_ot_env ( ispin ) % ener_x ( : ) & 01336 + qs_ot_env ( 1 ) % c_diis ( i ) * qs_ot_env ( ispin ) % ener_h_x ( i, : ) 01337 END DO 01338 ENDDO 01339 END IF 01340 qs_ot_env(1)%diis_iter=qs_ot_env(1)%diis_iter+1 01341 IF (qs_ot_env(1)%settings%safer_diis) THEN 01342 ! now, final check, is the step in fact in the direction of the -gradient ? 01343 ! if not we're walking towards a sadle point, and should avoid that 01344 ! the direction of the step is x_new-x_old 01345 tr_xold_gx=0.0_dp 01346 tr_xnew_gx=0.0_dp 01347 IF ( do_ks ) THEN 01348 DO ispin=1,nspin 01349 CALL cp_dbcsr_trace(qs_ot_env(ispin)%matrix_h_x(j)%matrix, & 01350 qs_ot_env(ispin)%matrix_gx, tmp,error=error) 01351 tr_xold_gx=tr_xold_gx+tmp 01352 CALL cp_dbcsr_trace(qs_ot_env(ispin)%matrix_x, & 01353 qs_ot_env(ispin)%matrix_gx, tmp,error=error) 01354 tr_xnew_gx=tr_xnew_gx+tmp 01355 IF (qs_ot_env(ispin)%settings%do_rotation) THEN 01356 CALL cp_dbcsr_trace(qs_ot_env(ispin)%rot_mat_h_x(j)%matrix, & 01357 qs_ot_env(ispin)%rot_mat_gx, tmp,error=error) 01358 tr_xold_gx=tr_xold_gx+0.5_dp*tmp 01359 CALL cp_dbcsr_trace(qs_ot_env(ispin)%rot_mat_x, & 01360 qs_ot_env(ispin)%rot_mat_gx, tmp,error=error) 01361 tr_xnew_gx=tr_xnew_gx+0.5_dp*tmp 01362 ENDIF 01363 ENDDO 01364 END IF 01365 ! ***SCP 01366 IF (do_scp_dft) THEN 01367 tmp = DOT_PRODUCT ( qs_ot_env ( 1 ) % h_x ( j, : ), qs_ot_env ( 1 ) % gx ( : ) ) 01368 CALL mp_sum ( tmp, qs_ot_env ( 1 ) % scp_para_env % group ) 01369 tr_xold_gx = tr_xold_gx+tmp 01370 tmp = DOT_PRODUCT ( qs_ot_env ( 1 ) % x ( : ), qs_ot_env ( 1 ) % gx ( : ) ) 01371 CALL mp_sum ( tmp, qs_ot_env ( 1 ) % scp_para_env % group ) 01372 tr_xnew_gx = tr_xnew_gx+tmp 01373 END IF 01374 IF (do_scp_nddo) THEN 01375 CALL cp_dbcsr_trace ( qs_ot_env ( 1 ) % hx_mat ( j ) % matrix, qs_ot_env ( 1 ) % gxmat, tmp, & 01376 local_sum=.TRUE., error=error ) 01377 tr_xold_gx = tr_xold_gx+tmp 01378 CALL cp_dbcsr_trace ( qs_ot_env ( 1 ) % xmat, qs_ot_env ( 1 ) % gxmat, tmp, local_sum=.TRUE.,& 01379 error=error ) 01380 tr_xnew_gx = tr_xnew_gx+tmp 01381 END IF 01382 ! ***SCP 01383 IF (do_ener) THEN 01384 DO ispin=1,nspin 01385 tmp = DOT_PRODUCT ( qs_ot_env ( ispin ) % ener_h_x ( j, : ), qs_ot_env ( ispin ) % ener_gx ( : ) ) 01386 tr_xold_gx = tr_xold_gx+tmp 01387 tmp = DOT_PRODUCT ( qs_ot_env ( ispin ) % ener_x ( : ), qs_ot_env ( ispin ) % ener_gx ( : ) ) 01388 tr_xnew_gx = tr_xnew_gx+tmp 01389 ENDDO 01390 END IF 01391 overlap=(tr_xnew_gx-tr_xold_gx) 01392 ! OK, bad luck, take a SD step along the preconditioned gradient 01393 IF (overlap.GT.0.0_dp) THEN 01394 qs_ot_env(1)%OT_METHOD_FULL="OT SD" 01395 IF ( do_ks ) THEN 01396 DO ispin=1,nspin 01397 CALL cp_dbcsr_set(qs_ot_env(ispin)%matrix_x,0.0_dp,error=error) 01398 CALL cp_dbcsr_add(qs_ot_env(ispin)%matrix_x, & 01399 qs_ot_env(ispin)%matrix_h_e(j)%matrix,& 01400 1.0_dp,1.0_dp,error=error) 01401 CALL cp_dbcsr_add(qs_ot_env(ispin)%matrix_x, & 01402 qs_ot_env(ispin)%matrix_h_x(j)%matrix,& 01403 1.0_dp,1.0_dp,error=error) 01404 IF (qs_ot_env(ispin)%settings%do_rotation) THEN 01405 CALL cp_dbcsr_set(qs_ot_env(ispin)%rot_mat_x,0.0_dp,error=error) 01406 CALL cp_dbcsr_add(qs_ot_env(ispin)%rot_mat_x, & 01407 qs_ot_env(ispin)%rot_mat_h_e(j)%matrix,& 01408 1.0_dp,1.0_dp,error=error) 01409 CALL cp_dbcsr_add(qs_ot_env(ispin)%rot_mat_x, & 01410 qs_ot_env(ispin)%rot_mat_h_x(j)%matrix,& 01411 1.0_dp,1.0_dp,error=error) 01412 ENDIF 01413 ENDDO 01414 END IF 01415 ! ***SCP 01416 IF (do_scp_dft) THEN 01417 qs_ot_env ( 1 ) % x ( : ) = 0._dp 01418 qs_ot_env ( 1 ) % x ( : ) = qs_ot_env ( 1 ) % x ( : ) + qs_ot_env ( 1 ) % h_e ( j, : ) 01419 qs_ot_env ( 1 ) % x ( : ) = qs_ot_env ( 1 ) % x ( : ) + qs_ot_env ( 1 ) % h_x ( j, : ) 01420 END IF 01421 IF (do_scp_nddo) THEN 01422 CALL cp_dbcsr_set ( qs_ot_env ( 1 ) % xmat , 0._dp , error=error) 01423 CALL cp_dbcsr_add ( qs_ot_env ( 1 ) % xmat, & 01424 qs_ot_env ( 1 ) % he_mat ( j ) % matrix, & 01425 alpha_scalar=1.0_dp, beta_scalar=1.0_dp, error=error ) 01426 CALL cp_dbcsr_add ( qs_ot_env ( 1 ) % xmat, qs_ot_env ( 1 ) % hx_mat ( j ) % matrix, & 01427 alpha_scalar=1.0_dp, beta_scalar=1.0_dp, error=error ) 01428 END IF 01429 ! ***SCP 01430 IF (do_ener) THEN 01431 DO ispin=1,nspin 01432 qs_ot_env ( ispin ) % ener_x ( : ) = 0._dp 01433 qs_ot_env ( ispin ) % ener_x ( : ) = qs_ot_env ( ispin ) % ener_x ( : ) + qs_ot_env ( ispin ) % ener_h_e ( j, : ) 01434 qs_ot_env ( ispin ) % ener_x ( : ) = qs_ot_env ( ispin ) % ener_x ( : ) + qs_ot_env ( ispin ) % ener_h_x ( j, : ) 01435 ENDDO 01436 END IF 01437 ENDIF 01438 ENDIF 01439 END SUBROUTINE ot_diis_step 01440 01441 ! ***************************************************************************** 01447 SUBROUTINE ot_broyden_step(qs_ot_env,error) 01448 TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env 01449 TYPE(cp_error_type), INTENT(inout) :: error 01450 01451 CHARACTER(len=*), PARAMETER :: routineN = 'ot_broyden_step', 01452 routineP = moduleN//':'//routineN 01453 01454 INTEGER :: diis_bound, diis_m, i, ispin, 01455 itmp, j, k, n, nspin 01456 INTEGER, ALLOCATABLE, DIMENSION(:) :: circ_index 01457 LOGICAL :: adaptive_sigma, do_ener, 01458 do_ks, do_scp_dft, 01459 do_scp_nddo, enable_flip, 01460 forget_history 01461 REAL(KIND=dp) :: beta, eta, gamma, omega, 01462 sigma, sigma_dec, sigma_min, 01463 tmp, tmp2 01464 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f, x 01465 REAL(KIND=dp), ALLOCATABLE, 01466 DIMENSION(:, :) :: G, S 01467 TYPE(cp_logger_type), POINTER :: logger 01468 TYPE(dbcsr_error_type) :: dbcsr_error 01469 01470 eta = qs_ot_env ( 1 ) % settings % broyden_eta 01471 omega = qs_ot_env ( 1 ) % settings % broyden_omega 01472 sigma_dec = qs_ot_env ( 1 ) % settings % broyden_sigma_decrease 01473 sigma_min = qs_ot_env ( 1 ) % settings % broyden_sigma_min 01474 forget_history = qs_ot_env ( 1 ) % settings % broyden_forget_history 01475 adaptive_sigma = qs_ot_env ( 1 ) % settings % broyden_adaptive_sigma 01476 enable_flip = qs_ot_env ( 1 ) % settings % broyden_enable_flip 01477 01478 do_ks = qs_ot_env ( 1 ) % settings % ks 01479 do_scp_dft = qs_ot_env ( 1 ) % settings % scp_dft 01480 do_scp_nddo = qs_ot_env ( 1 ) % settings % scp_nddo 01481 do_ener = qs_ot_env ( 1 ) % settings % do_ener 01482 01483 beta=qs_ot_env ( 1 ) % settings % broyden_beta 01484 gamma=qs_ot_env ( 1 ) % settings % broyden_gamma 01485 IF (adaptive_sigma) THEN 01486 IF (qs_ot_env ( 1 ) % broyden_adaptive_sigma .LT. 0.0_dp) THEN 01487 sigma=qs_ot_env ( 1 ) % settings % broyden_sigma 01488 ELSE 01489 sigma=qs_ot_env ( 1 ) % broyden_adaptive_sigma 01490 ENDIF 01491 ELSE 01492 sigma=qs_ot_env ( 1 ) % settings % broyden_sigma 01493 ENDIF 01494 01495 ! simplify our life.... 01496 IF (do_ener .OR. do_scp_nddo .OR. do_scp_dft .OR. .NOT. do_ks .OR. & 01497 qs_ot_env(1)%settings%do_rotation) THEN 01498 CALL stop_program(routineN,moduleN,__LINE__,"Not yet implemented") 01499 ENDIF 01500 ! 01501 nspin=SIZE(qs_ot_env) 01502 01503 diis_m=qs_ot_env(1)%settings%diis_m 01504 01505 IF (qs_ot_env(1)%diis_iter.lt.diis_m) THEN 01506 diis_bound=qs_ot_env(1)%diis_iter+1 01507 ELSE 01508 diis_bound=diis_m 01509 ENDIF 01510 01511 ! We want x:s, f:s and one random vector 01512 k = 2*diis_bound+1 01513 ALLOCATE(S(k,k)) 01514 ALLOCATE(G(k,k)) 01515 ALLOCATE(f(k)) 01516 ALLOCATE(x(k)) 01517 ALLOCATE(circ_index(diis_bound)) 01518 G = 0.0 01519 DO i=1,k 01520 G(i,i)=sigma 01521 ENDDO 01522 S = 0.0 01523 01524 01525 j = MOD(qs_ot_env(1)%diis_iter,diis_m)+1 ! index in the circular array 01526 01527 DO ispin=1,nspin 01528 CALL cp_dbcsr_copy(qs_ot_env(ispin)%matrix_h_x(j)%matrix,qs_ot_env(ispin)%matrix_x,error=error) 01529 ENDDO 01530 01531 IF (ASSOCIATED(qs_ot_env(1)%preconditioner)) THEN 01532 qs_ot_env(1)%gnorm=0.0_dp 01533 DO ispin=1,nspin 01534 CALL apply_preconditioner(qs_ot_env(ispin)%preconditioner, & 01535 qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_h_e(j)%matrix,error=error) 01536 CALL cp_dbcsr_trace(qs_ot_env(ispin)%matrix_gx,qs_ot_env(ispin)%matrix_h_e(j)%matrix, & 01537 tmp,error=error) 01538 qs_ot_env(1)%gnorm=qs_ot_env(1)%gnorm+tmp 01539 ENDDO 01540 IF (qs_ot_env(1)%gnorm .LT. 0.0_dp) THEN 01541 WRITE(cp_logger_get_default_unit_nr(logger),*) "WARNING Preconditioner not positive definite !" 01542 ENDIF 01543 DO ispin=1,nspin 01544 CALL cp_dbcsr_scale(qs_ot_env(ispin)%matrix_h_e(j)%matrix,-1.0_dp,error=error) 01545 ENDDO 01546 ELSE 01547 qs_ot_env(1)%gnorm=0.0_dp 01548 DO ispin=1,nspin 01549 CALL cp_dbcsr_trace(qs_ot_env(ispin)%matrix_gx,qs_ot_env(ispin)%matrix_gx,tmp,error=error) 01550 qs_ot_env(1)%gnorm=qs_ot_env(1)%gnorm+tmp 01551 CALL cp_dbcsr_add(qs_ot_env(ispin)%matrix_h_e(j)%matrix, & 01552 qs_ot_env(ispin)%matrix_gx,alpha_scalar=0.0_dp,beta_scalar=-1.0_dp,error=error) 01553 ENDDO 01554 ENDIF 01555 01556 k = 0 01557 n = 0 01558 CALL cp_dbcsr_get_info(qs_ot_env(1)%matrix_x,nfullrows_total=n) 01559 DO ispin=1,nspin 01560 CALL cp_dbcsr_get_info(qs_ot_env(ispin)%matrix_x,nfullcols_total=itmp) 01561 k=k+itmp 01562 ENDDO 01563 01564 ! Handling the case of no free variables to optimize 01565 IF (INT(n,KIND=int_8)*INT(k,KIND=int_8) /= 0) THEN 01566 qs_ot_env(1)%delta=SQRT(ABS(qs_ot_env(1)%gnorm)/(INT(n,KIND=int_8)*INT(k,KIND=int_8))) 01567 qs_ot_env(1)%gradient = - qs_ot_env(1)%gnorm 01568 ELSE 01569 qs_ot_env(1)%delta = 0.0_dp 01570 qs_ot_env(1)%gradient = 0.0_dp 01571 END IF 01572 01573 IF (diis_bound == diis_m) THEN 01574 DO i=1,diis_bound 01575 circ_index(i) = MOD(j+i-1,diis_m)+1 01576 ENDDO 01577 ELSE 01578 DO i=1,diis_bound 01579 circ_index(i) = i 01580 ENDDO 01581 ENDIF 01582 01583 S = 0.0_dp 01584 DO ispin=1,nspin 01585 CALL dbcsr_init_random(qs_ot_env(ispin)%matrix_x%matrix, dbcsr_error) 01586 DO i=1,diis_bound 01587 CALL cp_dbcsr_trace( & 01588 qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix, & 01589 qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix,tmp,error=error) 01590 S(i,i) = S(i,i) + tmp 01591 CALL cp_dbcsr_trace( & 01592 qs_ot_env(ispin)%matrix_h_e(circ_index(i))%matrix, & 01593 qs_ot_env(ispin)%matrix_h_e(circ_index(i))%matrix,tmp,error=error) 01594 S(i+diis_bound,i+diis_bound) = S(i+diis_bound,i+diis_bound) +tmp 01595 CALL cp_dbcsr_trace( & 01596 qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix, & 01597 qs_ot_env(ispin)%matrix_x,tmp,error=error) 01598 S(i,2*diis_bound+1) = S(i,2*diis_bound+1) + tmp 01599 S(i,2*diis_bound+1) = S(2*diis_bound+1,i) 01600 CALL cp_dbcsr_trace( & 01601 qs_ot_env(ispin)%matrix_h_e(circ_index(i))%matrix, & 01602 qs_ot_env(ispin)%matrix_x,tmp,error=error) 01603 S(i+diis_bound,2*diis_bound+1) = S(i+diis_bound,2*diis_bound+1) + tmp 01604 S(i+diis_bound,2*diis_bound+1) = S(2*diis_bound+1,diis_bound+i) 01605 DO k=(i+1),diis_bound 01606 CALL cp_dbcsr_trace( & 01607 qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix, & 01608 qs_ot_env(ispin)%matrix_h_x(circ_index(k))%matrix, & 01609 tmp,error=error) 01610 S(i,k) = S(i,k) + tmp 01611 S(k,i) = S(i,k) 01612 CALL cp_dbcsr_trace( & 01613 qs_ot_env(ispin)%matrix_h_e(circ_index(i))%matrix, & 01614 qs_ot_env(ispin)%matrix_h_e(circ_index(k))%matrix, & 01615 tmp,error=error) 01616 S(diis_bound+i,diis_bound+k) = S(diis_bound+i,diis_bound+k) + tmp 01617 S(diis_bound+k,diis_bound+i) = S(diis_bound+i,diis_bound+k) 01618 ENDDO 01619 DO k=1,diis_bound 01620 CALL cp_dbcsr_trace( & 01621 qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix, & 01622 qs_ot_env(ispin)%matrix_h_e(circ_index(k))%matrix,tmp,error=error) 01623 S(i,k+diis_bound) = S(i,k+diis_bound) + tmp 01624 S(k+diis_bound,i) = S(i,k+diis_bound) 01625 ENDDO 01626 ENDDO 01627 CALL cp_dbcsr_trace(qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_x,tmp,error=error) 01628 S(2*diis_bound+1,2*diis_bound+1) = S(2*diis_bound+1,2*diis_bound+1) + tmp 01629 ENDDO 01630 01631 ! normalize 01632 k = 2*diis_bound+1 01633 tmp = SQRT(S(k,k)) 01634 S(k,:) = S(k,:)/tmp 01635 S(:,k) = S(:,k)/tmp 01636 01637 IF (diis_bound .GT. 1) THEN 01638 tmp = 0.0_dp 01639 tmp2 = 0.0_dp 01640 i = diis_bound 01641 DO ispin=1,nspin 01642 ! dot product of differences 01643 CALL cp_dbcsr_trace( & 01644 qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix, & 01645 qs_ot_env(ispin)%matrix_h_e(circ_index(i))%matrix, & 01646 tmp,error=error) 01647 tmp2 = tmp2+tmp 01648 CALL cp_dbcsr_trace( & 01649 qs_ot_env(ispin)%matrix_h_x(circ_index(i-1))%matrix, & 01650 qs_ot_env(ispin)%matrix_h_e(circ_index(i))%matrix, & 01651 tmp,error=error) 01652 tmp2 = tmp2-tmp 01653 CALL cp_dbcsr_trace( & 01654 qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix, & 01655 qs_ot_env(ispin)%matrix_h_e(circ_index(i-1))%matrix, & 01656 tmp,error=error) 01657 tmp2 = tmp2-tmp 01658 CALL cp_dbcsr_trace( & 01659 qs_ot_env(ispin)%matrix_h_x(circ_index(i-1))%matrix, & 01660 qs_ot_env(ispin)%matrix_h_e(circ_index(i-1))%matrix, & 01661 tmp,error=error) 01662 tmp2 = tmp2+tmp 01663 ENDDO 01664 qs_ot_env(1)%c_broy(i-1) = tmp2 01665 ENDIF 01666 01667 01668 qs_ot_env(1)%energy_h(j) = qs_ot_env(1)%etotal 01669 01670 ! If we went uphill, do backtracking line search 01671 i = MINLOC(qs_ot_env(1)%energy_h(1:diis_bound),dim=1) 01672 IF (i .NE. j) THEN 01673 sigma = sigma_dec * sigma 01674 qs_ot_env(1)%OT_METHOD_FULL="OT BTRK" 01675 DO ispin=1,nspin 01676 CALL cp_dbcsr_set(qs_ot_env(ispin)%matrix_x,0.0_dp,error=error) 01677 CALL cp_dbcsr_add(qs_ot_env(ispin)%matrix_x, & 01678 qs_ot_env(ispin)%matrix_h_x(i)%matrix,& 01679 alpha_scalar=1.0_dp,beta_scalar=(1.0_dp-gamma),error=error) 01680 CALL cp_dbcsr_add(qs_ot_env(ispin)%matrix_x, & 01681 qs_ot_env(ispin)%matrix_h_x(circ_index(diis_bound))%matrix,& 01682 alpha_scalar=1.0_dp,beta_scalar=gamma,error=error) 01683 ENDDO 01684 ELSE 01685 ! Construct G 01686 DO i=2,diis_bound 01687 f = 0.0 01688 x = 0.0 01689 ! f is df_i 01690 x(i) = 1.0 01691 x(i-1) = -1.0 01692 ! x is dx_i 01693 f(diis_bound+i) = 1.0 01694 f(diis_bound+i-1) = -1.0 01695 tmp = 1.0_dp 01696 ! We want a pos def Hessian 01697 IF (enable_flip) THEN 01698 IF(qs_ot_env(1)%c_broy(i-1) .GT. 0) THEN 01699 !qs_ot_env(1)%OT_METHOD_FULL="OT FLIP" 01700 tmp = -1.0_dp 01701 ENDIF 01702 ENDIF 01703 01704 ! get dx-Gdf 01705 x = tmp*x - MATMUL(G, f) 01706 ! dfSdf 01707 ! we calculate matmul(S, f) twice. They're small... 01708 tmp = DOT_PRODUCT(f, MATMUL(S, f)) 01709 ! NOTE THAT S IS SYMMETRIC !!! 01710 f = MATMUL(S, f)/tmp 01711 ! the spread is an outer vector product 01712 G = G + SPREAD(x,dim=2,ncopies=SIZE(f))*SPREAD(f,dim=1,ncopies=SIZE(x)) 01713 ENDDO 01714 f = 0.0_dp 01715 f(2*diis_bound)=1.0_dp 01716 x = -beta*MATMUL(G,f) 01717 01718 ! OK, add the vectors now, this sums up to the proposed step 01719 DO ispin=1,nspin 01720 CALL cp_dbcsr_set(qs_ot_env(ispin)%matrix_x,0.0_dp,error=error) 01721 DO i=1, diis_bound 01722 CALL cp_dbcsr_add(qs_ot_env(ispin)%matrix_x, & 01723 qs_ot_env(ispin)%matrix_h_e(circ_index(i))%matrix,& 01724 alpha_scalar=1.0_dp,beta_scalar=-x(i+diis_bound),error=error) 01725 ENDDO 01726 DO i=1, diis_bound 01727 CALL cp_dbcsr_add(qs_ot_env(ispin)%matrix_x, & 01728 qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix,& 01729 alpha_scalar=1.0_dp,beta_scalar=x(i),error=error) 01730 ENDDO 01731 ENDDO 01732 01733 IF (adaptive_sigma) THEN 01734 tmp = new_sigma(G, S, diis_bound) 01735 !tmp = tmp * qs_ot_env ( 1 ) % settings % broyden_sigma 01736 tmp = tmp * eta 01737 sigma = MIN(omega * sigma, tmp) 01738 ENDIF 01739 01740 ! compute the inner product of direction of the step and gradient 01741 tmp = 0.0_dp 01742 DO ispin=1,nspin 01743 CALL cp_dbcsr_trace(qs_ot_env(ispin)%matrix_gx, & 01744 qs_ot_env(ispin)%matrix_x, & 01745 tmp2,error=error) 01746 tmp = tmp+tmp2 01747 ENDDO 01748 01749 DO ispin=1,nspin 01750 ! if the direction of the step is not in direction of the gradient, 01751 ! change step sign 01752 IF (tmp .GE. 0.0_dp) THEN 01753 qs_ot_env(1)%OT_METHOD_FULL="OT TURN" 01754 IF (forget_history) THEN 01755 qs_ot_env(1)%diis_iter=0 01756 ENDIF 01757 sigma = sigma*sigma_dec 01758 CALL cp_dbcsr_add(qs_ot_env(ispin)%matrix_x, & 01759 qs_ot_env(ispin)%matrix_h_x(circ_index(diis_bound))%matrix,& 01760 alpha_scalar=-1.0_dp,beta_scalar=1.0_dp,error=error) 01761 ELSE 01762 CALL cp_dbcsr_add(qs_ot_env(ispin)%matrix_x, & 01763 qs_ot_env(ispin)%matrix_h_x(circ_index(diis_bound))%matrix,& 01764 alpha_scalar=1.0_dp,beta_scalar=1.0_dp,error=error) 01765 ENDIF 01766 ENDDO 01767 ENDIF 01768 01769 01770 ! get rid of S, G, f, x, circ_index for next round 01771 DEALLOCATE(S,G,f,x,circ_index) 01772 01773 ! update for next round 01774 qs_ot_env(1)%diis_iter=qs_ot_env(1)%diis_iter+1 01775 qs_ot_env(1)%broyden_adaptive_sigma=MAX(sigma, sigma_min) 01776 01777 END SUBROUTINE ot_broyden_step 01778 01779 01780 FUNCTION new_sigma(G, S, n) RESULT(sigma) 01781 ! 01782 ! Calculate new sigma from eigenvalues of full size G by Arnoldi. 01783 ! 01784 ! ***************************************************************************** 01785 01786 01787 REAL(KIND=dp), DIMENSION(:, :) :: G, S 01788 INTEGER :: n 01789 REAL(KIND=dp) :: sigma 01790 01791 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigv 01792 REAL(KIND=dp), ALLOCATABLE, 01793 DIMENSION(:, :) :: H 01794 TYPE(cp_error_type) :: error 01795 01796 ALLOCATE(H(n,n)) 01797 CALL hess_G(G, S, H, n) 01798 ALLOCATE(eigv(n)) 01799 CALL diamat_all(H(1:n,1:n), eigv, error=error) 01800 01801 SELECT CASE(1) 01802 CASE(1) 01803 ! This estimator seems to work well. No theory. 01804 sigma = SUM(ABS(eigv**2))/SUM(ABS(eigv)) 01805 CASE(2) 01806 ! Estimator based on Frobenius norm minimizer 01807 sigma = SUM(ABS(eigv))/MAX(1, SIZE(eigv)) 01808 CASE(3) 01809 ! Estimator based on induced 2-norm 01810 sigma = (MAXVAL(ABS(eigv)) + MINVAL(ABS(eigv)))*0.5_dp 01811 END SELECT 01812 01813 DEALLOCATE(H, eigv) 01814 END FUNCTION new_sigma 01815 01816 SUBROUTINE hess_G(G, S, H, n) 01817 ! 01818 ! Make a hessenberg out of G into H. Cf Arnoldi. 01819 ! Inner product is weighted by S. 01820 ! Possible lucky breakdown at n. 01821 ! 01822 ! ***************************************************************************** 01823 REAL(KIND=dp), DIMENSION(:, :) :: G, S, H 01824 INTEGER :: n 01825 01826 INTEGER :: i, j, k 01827 REAL(KIND=dp) :: tmp 01828 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: v 01829 REAL(KIND=dp), ALLOCATABLE, 01830 DIMENSION(:, :) :: Q 01831 01832 i = SIZE(G, 1) 01833 k = SIZE(H, 1) 01834 ALLOCATE(Q(i,k)) 01835 ALLOCATE(v(i)) 01836 H = 0.0_dp 01837 Q = 0.0_dp 01838 01839 Q(:,1) = 1.0_dp 01840 tmp = SQRT(DOT_PRODUCT(Q(:,1), MATMUL(S, Q(:,1)))) 01841 Q = Q/tmp 01842 01843 01844 DO i=1,k 01845 v = MATMUL(G, Q(:,i)) 01846 DO j=1,i 01847 H(j,i) = DOT_PRODUCT(Q(:,j), MATMUL(S,v) ) 01848 v = v - H(j,i) * Q(:,j) 01849 ENDDO 01850 IF (i .LT. k) THEN 01851 tmp = DOT_PRODUCT(v, MATMUL(S, v)) 01852 IF (tmp .LE. 0.0_dp) THEN 01853 n = i 01854 EXIT 01855 ENDIF 01856 tmp = SQRT(tmp) 01857 ! Lucky breakdown 01858 IF (ABS(tmp) .LT. 1e-9_dp) THEN 01859 n = i 01860 EXIT 01861 ENDIF 01862 H(i+1,i) = tmp 01863 Q(:,i+1) = v/H(i+1,i) 01864 ENDIF 01865 ENDDO 01866 01867 DEALLOCATE(Q, v) 01868 END SUBROUTINE hess_G 01869 01870 END MODULE qs_ot_minimizer
1.7.3