1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief density matrix optimization using exponential transformations 8!> \par History 9!> 2012.05 created [Florian Schiffmann] 10!> \author Florian Schiffmann 11! ************************************************************************************************** 12 13MODULE dm_ls_scf_curvy 14 USE bibliography, ONLY: Shao2003,& 15 cite_reference 16 USE cp_log_handling, ONLY: cp_get_default_logger,& 17 cp_logger_get_default_unit_nr,& 18 cp_logger_type 19 USE dbcsr_api, ONLY: & 20 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_dot, dbcsr_filter, dbcsr_frobenius_norm, & 21 dbcsr_multiply, dbcsr_norm, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_transposed, & 22 dbcsr_type, dbcsr_type_no_symmetry 23 USE dm_ls_scf_types, ONLY: ls_scf_curvy_type,& 24 ls_scf_env_type 25 USE input_constants, ONLY: ls_scf_line_search_3point,& 26 ls_scf_line_search_3point_2d 27 USE iterate_matrix, ONLY: purify_mcweeny 28 USE kinds, ONLY: dp 29 USE machine, ONLY: m_flush 30 USE mathconstants, ONLY: ifac 31 USE mathlib, ONLY: invmat 32#include "./base/base_uses.f90" 33 34 IMPLICIT NONE 35 36 PRIVATE 37 38 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dm_ls_scf_curvy' 39 40 PUBLIC :: dm_ls_curvy_optimization, deallocate_curvy_data 41 42CONTAINS 43 44! ************************************************************************************************** 45!> \brief driver routine for Head-Gordon curvy step approach 46!> \param ls_scf_env ... 47!> \param energy ... 48!> \param check_conv ... 49!> \par History 50!> 2012.05 created [Florian Schiffmann] 51!> \author Florian Schiffmann 52! ************************************************************************************************** 53 54 SUBROUTINE dm_ls_curvy_optimization(ls_scf_env, energy, check_conv) 55 TYPE(ls_scf_env_type) :: ls_scf_env 56 REAL(KIND=dp) :: energy 57 LOGICAL :: check_conv 58 59 CHARACTER(LEN=*), PARAMETER :: routineN = 'dm_ls_curvy_optimization', & 60 routineP = moduleN//':'//routineN 61 62 INTEGER :: handle, i, lsstep 63 64 CALL timeset(routineN, handle) 65 66 CALL cite_reference(Shao2003) 67 68! Upon first call initialize all matrices needed curing optimization 69! In addtion transform P into orthonormal basis. Will be scaled by 0.5 in closed shell case 70! Only to be done once as it will be stored and reused afterwards 71! TRS4 might yield a non-idempotent P therefore McWeeny purification is applied on initial P 72 73 IF (.NOT. ALLOCATED(ls_scf_env%curvy_data%matrix_dp)) THEN 74 CALL init_curvy(ls_scf_env%curvy_data, ls_scf_env%matrix_s, ls_scf_env%nspins) 75 ls_scf_env%curvy_data%line_search_step = 1 76 77 IF (ls_scf_env%curvy_data%line_search_type == ls_scf_line_search_3point_2d) THEN 78 DO i = 1, ls_scf_env%nspins 79 CALL dbcsr_copy(ls_scf_env%curvy_data%matrix_psave(i, 1), & 80 ls_scf_env%matrix_p(i)) 81 END DO 82 END IF 83 IF (ls_scf_env%nspins == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(1), 0.5_dp) 84 CALL transform_matrix_orth(ls_scf_env%matrix_p, ls_scf_env%matrix_s_sqrt, & 85 ls_scf_env%eps_filter) 86 CALL purify_mcweeny(ls_scf_env%matrix_p, ls_scf_env%eps_filter, 3) 87 DO i = 1, ls_scf_env%nspins 88 CALL dbcsr_copy(ls_scf_env%curvy_data%matrix_p(i), ls_scf_env%matrix_p(i)) 89 END DO 90 END IF 91 92 lsstep = ls_scf_env%curvy_data%line_search_step 93 94! If new search direction has to be computed transform H into the orthnormal basis 95 96 IF (ls_scf_env%curvy_data%line_search_step == 1) & 97 CALL transform_matrix_orth(ls_scf_env%matrix_ks, ls_scf_env%matrix_s_sqrt_inv, & 98 ls_scf_env%eps_filter) 99 100! Set the energies for the line search and make sure to give the correct energy back to scf_main 101 ls_scf_env%curvy_data%energies(lsstep) = energy 102 IF (lsstep .NE. 1) energy = ls_scf_env%curvy_data%energies(1) 103 104! start the optimization by calling the driver routine or simply combine saved P(2D line search) 105 IF (lsstep .LE. 2) THEN 106 CALL optimization_step(ls_scf_env%curvy_data, ls_scf_env) 107 ELSE IF (lsstep == ls_scf_env%curvy_data%line_search_type) THEN 108! line_search type has the value appropriate to the number of energy calculations needed 109 CALL optimization_step(ls_scf_env%curvy_data, ls_scf_env) 110 ELSE 111 CALL new_p_from_save(ls_scf_env%matrix_p, ls_scf_env%curvy_data%matrix_psave, lsstep, & 112 ls_scf_env%curvy_data%double_step_size) 113 ls_scf_env%curvy_data%line_search_step = ls_scf_env%curvy_data%line_search_step + 1 114 CALL timestop(handle) 115 RETURN 116 END IF 117 lsstep = ls_scf_env%curvy_data%line_search_step 118 119! transform new density matrix back into nonorthonormal basis (again scaling might apply) 120 121 CALL transform_matrix_orth(ls_scf_env%matrix_p, ls_scf_env%matrix_s_sqrt_inv, & 122 ls_scf_env%eps_filter) 123 IF (ls_scf_env%nspins == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(1), 2.0_dp) 124 125! P-matrices only need to be stored in case of 2D line search 126 IF (lsstep .LE. 3 .AND. ls_scf_env%curvy_data%line_search_type == ls_scf_line_search_3point_2d) THEN 127 DO i = 1, ls_scf_env%nspins 128 CALL dbcsr_copy(ls_scf_env%curvy_data%matrix_psave(i, lsstep), & 129 ls_scf_env%matrix_p(i)) 130 END DO 131 END IF 132 check_conv = lsstep == 1 133 134 CALL timestop(handle) 135 136 END SUBROUTINE dm_ls_curvy_optimization 137 138! ************************************************************************************************** 139!> \brief low level routine for Head-Gordons curvy step approach 140!> computes gradients, performs a cg and line search, 141!> and evaluates the BCH series to obtain the new P matrix 142!> \param curvy_data ... 143!> \param ls_scf_env ... 144!> \par History 145!> 2012.05 created [Florian Schiffmann] 146!> \author Florian Schiffmann 147! ************************************************************************************************** 148 149 SUBROUTINE optimization_step(curvy_data, ls_scf_env) 150 TYPE(ls_scf_curvy_type) :: curvy_data 151 TYPE(ls_scf_env_type) :: ls_scf_env 152 153 CHARACTER(LEN=*), PARAMETER :: routineN = 'optimization_step', & 154 routineP = moduleN//':'//routineN 155 156 INTEGER :: handle, ispin 157 REAL(KIND=dp) :: filter, step_size(2) 158 159! Upon first line search step compute new search direction and apply CG if required 160 161 CALL timeset(routineN, handle) 162 163 IF (curvy_data%line_search_step == 1) THEN 164 curvy_data%step_size = MAXVAL(curvy_data%step_size) 165 curvy_data%step_size = MIN(MAX(0.10_dp, 0.5_dp*ABS(curvy_data%step_size(1))), 0.5_dp) 166! Dynamic eps_filter for newton steps 167 filter = MAX(ls_scf_env%eps_filter*curvy_data%min_filter, & 168 ls_scf_env%eps_filter*curvy_data%filter_factor) 169 CALL compute_direction_newton(curvy_data%matrix_p, ls_scf_env%matrix_ks, & 170 curvy_data%matrix_dp, filter, curvy_data%fix_shift, curvy_data%shift, & 171 curvy_data%cg_numer, curvy_data%cg_denom, curvy_data%min_shift) 172 curvy_data%filter_factor = curvy_data%scale_filter*curvy_data%filter_factor 173 step_size = curvy_data%step_size 174 curvy_data%BCH_saved = 0 175 ELSE IF (curvy_data%line_search_step == 2) THEN 176 step_size = curvy_data%step_size 177 IF (curvy_data%energies(1) - curvy_data%energies(2) .GT. 0.0_dp) THEN 178 curvy_data%step_size = curvy_data%step_size*2.0_dp 179 curvy_data%double_step_size = .TRUE. 180 ELSE 181 curvy_data%step_size = curvy_data%step_size*0.5_dp 182 curvy_data%double_step_size = .FALSE. 183 END IF 184 step_size = curvy_data%step_size 185 ELSE IF (curvy_data%line_search_step == ls_scf_line_search_3point_2d) THEN 186 CALL line_search_2d(curvy_data%energies, curvy_data%step_size) 187 step_size = curvy_data%step_size 188 ELSE IF (curvy_data%line_search_step == ls_scf_line_search_3point) THEN 189 CALL line_search_3pnt(curvy_data%energies, curvy_data%step_size) 190 step_size = curvy_data%step_size 191 END IF 192 193 CALL update_p_exp(curvy_data%matrix_p, ls_scf_env%matrix_p, curvy_data%matrix_dp, & 194 curvy_data%matrix_BCH, ls_scf_env%eps_filter, step_size, curvy_data%BCH_saved, & 195 curvy_data%n_bch_hist) 196 197! line_search type has the value appropriate to the numeber of energy calculations needed 198 curvy_data%line_search_step = MOD(curvy_data%line_search_step, curvy_data%line_search_type) + 1 199 IF (curvy_data%line_search_step == 1) THEN 200 DO ispin = 1, SIZE(curvy_data%matrix_p) 201 CALL dbcsr_copy(curvy_data%matrix_p(ispin), ls_scf_env%matrix_p(ispin)) 202 END DO 203 END IF 204 CALL timestop(handle) 205 206 END SUBROUTINE optimization_step 207 208! ************************************************************************************************** 209!> \brief Perform a 6pnt-2D line search for spin polarized calculations. 210!> Fit a 2D parabolic function to 6 points 211!> \param energies ... 212!> \param step_size ... 213!> \par History 214!> 2012.05 created [Florian Schiffmann] 215!> \author Florian Schiffmann 216! ************************************************************************************************** 217 218 SUBROUTINE line_search_2d(energies, step_size) 219 REAL(KIND=dp) :: energies(6), step_size(2) 220 221 CHARACTER(LEN=*), PARAMETER :: routineN = 'line_search_2d', routineP = moduleN//':'//routineN 222 223 INTEGER :: info, unit_nr 224 REAL(KIND=dp) :: e_pred, param(6), s1, s1sq, s2, s2sq, & 225 sys_lin_eq(6, 6), tmp_e, v1, v2 226 TYPE(cp_logger_type), POINTER :: logger 227 228 logger => cp_get_default_logger() 229 IF (energies(1) - energies(2) .LT. 0._dp) THEN 230 tmp_e = energies(2); energies(2) = energies(3); energies(3) = tmp_e 231 step_size = step_size*2.0_dp 232 END IF 233 IF (logger%para_env%ionode) THEN 234 unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) 235 ELSE 236 unit_nr = -1 237 ENDIF 238 s1 = 0.5_dp*step_size(1); s2 = step_size(1); s1sq = s1**2; s2sq = s2**2 239 sys_lin_eq = 0.0_dp; sys_lin_eq(:, 6) = 1.0_dp 240 sys_lin_eq(2, 1) = s1sq; sys_lin_eq(2, 2) = s1sq; sys_lin_eq(2, 3) = s1sq; sys_lin_eq(2, 4) = s1; sys_lin_eq(2, 5) = s1 241 sys_lin_eq(3, 1) = s2sq; sys_lin_eq(3, 2) = s2sq; sys_lin_eq(3, 3) = s2sq; sys_lin_eq(3, 4) = s2; sys_lin_eq(3, 5) = s2 242 sys_lin_eq(4, 3) = s1sq; sys_lin_eq(4, 5) = s1 243 sys_lin_eq(5, 1) = s1sq; sys_lin_eq(5, 4) = s1 244 sys_lin_eq(6, 3) = s2sq; sys_lin_eq(6, 5) = s2 245 246 CALL invmat(sys_lin_eq, info) 247 param = MATMUL(sys_lin_eq, energies) 248 v1 = (param(2)*param(4))/(2.0_dp*param(1)) - param(5) 249 v2 = -(param(2)**2)/(2.0_dp*param(1)) + 2.0_dp*param(3) 250 step_size(2) = v1/v2 251 step_size(1) = (-param(2)*step_size(2) - param(4))/(2.0_dp*param(1)) 252 IF (step_size(1) .LT. 0.0_dp) step_size(1) = 1.0_dp 253 IF (step_size(2) .LT. 0.0_dp) step_size(2) = 1.0_dp 254! step_size(1)=MIN(step_size(1),2.0_dp) 255! step_size(2)=MIN(step_size(2),2.0_dp) 256 e_pred = param(1)*step_size(1)**2 + param(2)*step_size(1)*step_size(2) + & 257 param(3)*step_size(2)**2 + param(4)*step_size(1) + param(5)*step_size(2) + param(6) 258 IF (unit_nr .GT. 0) WRITE (unit_nr, "(t3,a,F10.5,F10.5,A,F20.9)") & 259 " Line Search: Step Size", step_size, " Predicted energy", e_pred 260 e_pred = param(1)*s1**2 + param(2)*s2*s1*0.0_dp + & 261 param(3)*s1**2*0.0_dp + param(4)*s1 + param(5)*s1*0.0_dp + param(6) 262 263 END SUBROUTINE line_search_2d 264 265! ************************************************************************************************** 266!> \brief Perform a 3pnt line search 267!> \param energies ... 268!> \param step_size ... 269!> \par History 270!> 2012.05 created [Florian Schiffmann] 271!> \author Florian Schiffmann 272! ************************************************************************************************** 273 274 SUBROUTINE line_search_3pnt(energies, step_size) 275 REAL(KIND=dp) :: energies(3), step_size(2) 276 277 CHARACTER(LEN=*), PARAMETER :: routineN = 'line_search_3pnt', & 278 routineP = moduleN//':'//routineN 279 280 INTEGER :: unit_nr 281 REAL(KIND=dp) :: a, b, c, e_pred, min_val, step1, tmp, & 282 tmp_e 283 TYPE(cp_logger_type), POINTER :: logger 284 285 logger => cp_get_default_logger() 286 IF (energies(1) - energies(2) .LT. 0._dp) THEN 287 tmp_e = energies(2); energies(2) = energies(3); energies(3) = tmp_e 288 step_size = step_size*2.0_dp 289 END IF 290 IF (logger%para_env%ionode) THEN 291 unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) 292 ELSE 293 unit_nr = -1 294 ENDIF 295 step1 = 0.5_dp*step_size(1) 296 c = energies(1) 297 a = (energies(3) + c - 2.0_dp*energies(2))/(2.0_dp*step1**2) 298 b = (energies(2) - c - a*step1**2)/step1 299 IF (a .LT. 1.0E-12_dp) a = -1.0E-12_dp 300 min_val = -b/(2.0_dp*a) 301 e_pred = a*min_val**2 + b*min_val + c 302 tmp = step_size(1) 303 IF (e_pred .LT. energies(1) .AND. e_pred .LT. energies(2)) THEN 304 step_size = MAX(-1.0_dp, & 305 MIN(min_val, 10_dp*step_size)) 306 ELSE 307 step_size = 1.0_dp 308 END IF 309 e_pred = a*(step_size(1))**2 + b*(step_size(1)) + c 310 IF (unit_nr .GT. 0) THEN 311 WRITE (unit_nr, "(t3,a,f16.8,a,F20.9)") "Line Search: Step Size", step_size(1), " Predicted energy", e_pred 312 CALL m_flush(unit_nr) 313 ENDIF 314 END SUBROUTINE line_search_3pnt 315 316! ************************************************************************************************** 317!> \brief Get a new search direction. Iterate to obtain a Newton like step 318!> Refine with a CG update of the search direction 319!> \param matrix_p ... 320!> \param matrix_ks ... 321!> \param matrix_dp ... 322!> \param eps_filter ... 323!> \param fix_shift ... 324!> \param curvy_shift ... 325!> \param cg_numer ... 326!> \param cg_denom ... 327!> \param min_shift ... 328!> \par History 329!> 2012.05 created [Florian Schiffmann] 330!> \author Florian Schiffmann 331! ************************************************************************************************** 332 333 SUBROUTINE compute_direction_newton(matrix_p, matrix_ks, matrix_dp, eps_filter, fix_shift, & 334 curvy_shift, cg_numer, cg_denom, min_shift) 335 TYPE(dbcsr_type), DIMENSION(:) :: matrix_p, matrix_ks, matrix_dp 336 REAL(KIND=dp) :: eps_filter 337 LOGICAL :: fix_shift(2) 338 REAL(KIND=dp) :: curvy_shift(2), cg_numer(2), & 339 cg_denom(2), min_shift 340 341 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_direction_newton', & 342 routineP = moduleN//':'//routineN 343 344 INTEGER :: handle, i, ispin, ncyc, nspin, unit_nr 345 LOGICAL :: at_limit 346 REAL(KIND=dp) :: beta, conv_val, maxel, old_conv, shift 347 TYPE(cp_logger_type), POINTER :: logger 348 TYPE(dbcsr_type) :: matrix_Ax, matrix_b, matrix_cg, & 349 matrix_dp_old, matrix_PKs, matrix_res, & 350 matrix_tmp, matrix_tmp1 351 352 logger => cp_get_default_logger() 353 354 IF (logger%para_env%ionode) THEN 355 unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) 356 ELSE 357 unit_nr = -1 358 ENDIF 359 CALL timeset(routineN, handle) 360 nspin = SIZE(matrix_p) 361 362 CALL dbcsr_create(matrix_PKs, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry) 363 CALL dbcsr_create(matrix_Ax, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry) 364 CALL dbcsr_create(matrix_tmp, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry) 365 CALL dbcsr_create(matrix_tmp1, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry) 366 CALL dbcsr_create(matrix_res, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry) 367 CALL dbcsr_create(matrix_cg, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry) 368 CALL dbcsr_create(matrix_b, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry) 369 CALL dbcsr_create(matrix_dp_old, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry) 370 371 DO ispin = 1, nspin 372 CALL dbcsr_copy(matrix_dp_old, matrix_dp(ispin)) 373 374! Precompute some matrices to save work during iterations 375 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p(ispin), matrix_ks(ispin), & 376 0.0_dp, matrix_PKs, filter_eps=eps_filter) 377 CALL dbcsr_transposed(matrix_b, matrix_PKs) 378 CALL dbcsr_copy(matrix_cg, matrix_b) 379 380! Starting CG with guess 0-matrix gives -2*gradient=[Ks*P-(Ks*P)T] for cg_matrix in second step 381 CALL dbcsr_add(matrix_cg, matrix_PKs, 2.0_dp, -2.0_dp) 382 383! Residual matrix in first step=cg matrix. Keep Pks for later use in CG! 384 CALL dbcsr_copy(matrix_res, matrix_cg) 385 386! Precompute -FP-[FP]T which will be used throughout the CG iterations 387 CALL dbcsr_add(matrix_b, matrix_PKs, -1.0_dp, -1.0_dp) 388 389! Setup some values to check convergence and safty checks for eigenvalue shifting 390 CALL dbcsr_norm(matrix_res, which_norm=2, norm_scalar=old_conv) 391 old_conv = dbcsr_frobenius_norm(matrix_res) 392 shift = MIN(10.0_dp, MAX(min_shift, 0.05_dp*old_conv)) 393 conv_val = MAX(0.010_dp*old_conv, 100.0_dp*eps_filter) 394 old_conv = 100.0_dp 395 IF (fix_shift(ispin)) THEN 396 shift = MAX(min_shift, MIN(10.0_dp, MAX(shift, curvy_shift(ispin) - 0.5_dp*curvy_shift(ispin)))) 397 curvy_shift(ispin) = shift 398 END IF 399 400! Begin the real optimization loop 401 CALL dbcsr_set(matrix_dp(ispin), 0.0_dp) 402 ncyc = 10 403 DO i = 1, ncyc 404 405! One step to compute: -FPD-DPF-DFP-PFD (not obvious but symmetry allows for some tricks) 406 CALL commutator_symm(matrix_b, matrix_cg, matrix_Ax, eps_filter, 1.0_dp) 407 408! Compute the missing bits 2*(FDP+PDF) (again use symmetry to compute as a commutator) 409 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_cg, matrix_p(ispin), & 410 0.0_dp, matrix_tmp, filter_eps=eps_filter) 411 CALL commutator_symm(matrix_ks(ispin), matrix_tmp, matrix_tmp1, eps_filter, 2.0_dp) 412 CALL dbcsr_add(matrix_Ax, matrix_tmp1, 1.0_dp, 1.0_dp) 413 414! Apply the shift and hope it's enough to stabilize the CG iterations 415 CALL dbcsr_add(matrix_Ax, matrix_cg, 1.0_dp, shift) 416 417 CALL compute_cg_matrices(matrix_Ax, matrix_res, matrix_cg, matrix_dp(ispin), & 418 matrix_tmp, eps_filter, at_limit) 419 CALL dbcsr_filter(matrix_cg, eps_filter) 420 421! check for convergence of the newton step 422 maxel = dbcsr_frobenius_norm(matrix_res) 423 IF (unit_nr .GT. 0) THEN 424 WRITE (unit_nr, "(T3,A,F12.6)") "Convergence of Newton iteration ", maxel 425 CALL m_flush(unit_nr) 426 ENDIF 427 at_limit = at_limit .OR. (old_conv/maxel .LT. 1.01_dp) 428 old_conv = maxel 429 IF (i == ncyc .AND. maxel/conv_val .GT. 5.0_dp) THEN 430 fix_shift(ispin) = .TRUE. 431 curvy_shift(ispin) = 4.0_dp*shift 432 END IF 433 IF (maxel .LT. conv_val .OR. at_limit) EXIT 434 END DO 435 436! Refine the Newton like search direction with a preconditioned cg update 437 CALL dbcsr_transposed(matrix_b, matrix_PKs) 438 !compute b= -2*KsP+2*PKs=-(2*gradient) 439 CALL dbcsr_copy(matrix_cg, matrix_b) 440 CALL dbcsr_add(matrix_cg, matrix_PKs, 1.0_dp, -1.0_dp) 441 cg_denom(ispin) = cg_numer(ispin) 442 CALL dbcsr_dot(matrix_cg, matrix_dp(ispin), cg_numer(ispin)) 443 beta = cg_numer(ispin)/MAX(cg_denom(ispin), 1.0E-6_dp) 444 IF (beta .LT. 1.0_dp) THEN 445 beta = MAX(0.0_dp, beta) 446 CALL dbcsr_add(matrix_dp(ispin), matrix_dp_old, 1.0_dp, beta) 447 END IF 448 IF (unit_nr .GT. 0) WRITE (unit_nr, "(A)") " " 449 END DO 450 451 CALL dbcsr_release(matrix_PKs) 452 CALL dbcsr_release(matrix_dp_old) 453 CALL dbcsr_release(matrix_b) 454 CALL dbcsr_release(matrix_Ax) 455 CALL dbcsr_release(matrix_tmp) 456 CALL dbcsr_release(matrix_tmp1) 457 CALL dbcsr_release(matrix_b) 458 CALL dbcsr_release(matrix_res) 459 CALL dbcsr_release(matrix_cg) 460 461 IF (unit_nr .GT. 0) CALL m_flush(unit_nr) 462 CALL timestop(handle) 463 END SUBROUTINE compute_direction_newton 464 465! ************************************************************************************************** 466!> \brief compute the optimal step size of the current cycle and update the 467!> matrices needed to solve the system of linear equations 468!> \param Ax ... 469!> \param res ... 470!> \param cg ... 471!> \param deltp ... 472!> \param tmp ... 473!> \param eps_filter ... 474!> \param at_limit ... 475!> \par History 476!> 2012.05 created [Florian Schiffmann] 477!> \author Florian Schiffmann 478! ************************************************************************************************** 479 480 SUBROUTINE compute_cg_matrices(Ax, res, cg, deltp, tmp, eps_filter, at_limit) 481 TYPE(dbcsr_type) :: Ax, res, cg, deltp, tmp 482 REAL(KIND=dp) :: eps_filter 483 LOGICAL :: at_limit 484 485 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_cg_matrices', & 486 routineP = moduleN//':'//routineN 487 488 INTEGER :: i, info 489 REAL(KIND=dp) :: alpha, beta, devi(3), fac, fac1, & 490 lin_eq(3, 3), new_norm, norm_cA, & 491 norm_rr, vec(3) 492 493 at_limit = .FALSE. 494 CALL dbcsr_dot(res, res, norm_rr) 495 CALL dbcsr_dot(cg, Ax, norm_cA) 496 lin_eq = 0.0_dp 497 fac = norm_rr/norm_cA 498 fac1 = fac 499! Use a 3point line serach and a fit to a quadratic function to determine optimal step size 500 DO i = 1, 3 501 CALL dbcsr_copy(tmp, res) 502 CALL dbcsr_add(tmp, Ax, 1.0_dp, -fac) 503 devi(i) = dbcsr_frobenius_norm(tmp) 504 lin_eq(i, :) = (/fac**2, fac, 1.0_dp/) 505 fac = fac1 + fac1*((-1)**i)*0.5_dp 506 END DO 507 CALL invmat(lin_eq, info) 508 vec = MATMUL(lin_eq, devi) 509 alpha = -vec(2)/(2.0_dp*vec(1)) 510 fac = SQRT(norm_rr/(norm_cA*alpha)) 511!scale the previous matrices to match the step size 512 CALL dbcsr_scale(Ax, fac) 513 CALL dbcsr_scale(cg, fac) 514 norm_cA = norm_cA*fac**2 515 516! USe CG to get the new matrices 517 alpha = norm_rr/norm_cA 518 CALL dbcsr_add(res, Ax, 1.0_dp, -alpha) 519 CALL dbcsr_dot(res, res, new_norm) 520 IF (norm_rr .LT. eps_filter*0.001_dp .OR. new_norm .LT. eps_filter*0.001_dp) THEN 521 beta = 0.0_dp 522 at_limit = .TRUE. 523 ELSE 524 beta = new_norm/norm_rr 525 CALL dbcsr_add(deltp, cg, 1.0_dp, alpha) 526 END IF 527 beta = new_norm/norm_rr 528 CALL dbcsr_add(cg, res, beta, 1.0_dp) 529 530 END SUBROUTINE compute_cg_matrices 531 532! ************************************************************************************************** 533!> \brief Only for 2D line serach. Use saved P-components to construct new 534!> test density matrix. Takes care as well, whether step_size 535!> increased or decreased during 2nd step and combines matrices accordingly 536!> \param matrix_p ... 537!> \param matrix_psave ... 538!> \param lsstep ... 539!> \param DOUBLE ... 540!> \par History 541!> 2012.05 created [Florian Schiffmann] 542!> \author Florian Schiffmann 543! ************************************************************************************************** 544 545 SUBROUTINE new_p_from_save(matrix_p, matrix_psave, lsstep, DOUBLE) 546 TYPE(dbcsr_type), DIMENSION(:) :: matrix_p 547 TYPE(dbcsr_type), DIMENSION(:, :) :: matrix_psave 548 INTEGER :: lsstep 549 LOGICAL :: DOUBLE 550 551 SELECT CASE (lsstep) 552 CASE (3) 553 CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 1)) 554 IF (DOUBLE) THEN 555 CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 2)) 556 ELSE 557 CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 3)) 558 END IF 559 CASE (4) 560 IF (DOUBLE) THEN 561 CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 2)) 562 ELSE 563 CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 3)) 564 END IF 565 CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 1)) 566 CASE (5) 567 CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 1)) 568 IF (DOUBLE) THEN 569 CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 3)) 570 ELSE 571 CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 2)) 572 END IF 573 END SELECT 574 575 END SUBROUTINE new_p_from_save 576 577! ************************************************************************************************** 578!> \brief computes a commutator exploiting symmetry RES=k*[A,B]=k*[AB-(AB)T] 579!> \param a ... 580!> \param b ... 581!> \param res ... 582!> \param eps_filter filtering threshold for sparse matrices 583!> \param prefac prefactor k in above equation 584!> \par History 585!> 2012.05 created [Florian Schiffmann] 586!> \author Florian Schiffmann 587! ************************************************************************************************** 588 589 SUBROUTINE commutator_symm(a, b, res, eps_filter, prefac) 590 TYPE(dbcsr_type) :: a, b, res 591 REAL(KIND=dp) :: eps_filter, prefac 592 593 CHARACTER(LEN=*), PARAMETER :: routineN = 'commutator_symm', & 594 routineP = moduleN//':'//routineN 595 596 INTEGER :: handle 597 TYPE(dbcsr_type) :: work 598 599 CALL timeset(routineN, handle) 600 601 CALL dbcsr_create(work, template=a, matrix_type=dbcsr_type_no_symmetry) 602 603 CALL dbcsr_multiply("N", "N", prefac, a, b, 0.0_dp, res, filter_eps=eps_filter) 604 CALL dbcsr_transposed(work, res) 605 CALL dbcsr_add(res, work, 1.0_dp, -1.0_dp) 606 607 CALL dbcsr_release(work) 608 609 CALL timestop(handle) 610 END SUBROUTINE commutator_symm 611 612! ************************************************************************************************** 613!> \brief Use the BCH update to get the new idempotent P 614!> Numerics don't allow for perfect idempotency, therefore a mc weeny 615!> step is used to make sure we stay close to the idempotent surface 616!> \param matrix_p_in ... 617!> \param matrix_p_out ... 618!> \param matrix_dp ... 619!> \param matrix_BCH ... 620!> \param threshold ... 621!> \param step_size ... 622!> \param BCH_saved ... 623!> \param n_bch_hist ... 624!> \par History 625!> 2012.05 created [Florian Schiffmann] 626!> \author Florian Schiffmann 627! ************************************************************************************************** 628 629 SUBROUTINE update_p_exp(matrix_p_in, matrix_p_out, matrix_dp, matrix_BCH, threshold, step_size, & 630 BCH_saved, n_bch_hist) 631 TYPE(dbcsr_type), DIMENSION(:) :: matrix_p_in, matrix_p_out, matrix_dp 632 TYPE(dbcsr_type), DIMENSION(:, :) :: matrix_BCH 633 REAL(KIND=dp) :: threshold, step_size(2) 634 INTEGER :: BCH_saved(2), n_bch_hist 635 636 CHARACTER(LEN=*), PARAMETER :: routineN = 'update_p_exp', routineP = moduleN//':'//routineN 637 638 INTEGER :: handle, i, ispin, nsave, nspin, unit_nr 639 LOGICAL :: save_BCH 640 REAL(KIND=dp) :: frob_norm, step_fac 641 TYPE(cp_logger_type), POINTER :: logger 642 TYPE(dbcsr_type) :: matrix, matrix_tmp 643 644 CALL timeset(routineN, handle) 645 646 logger => cp_get_default_logger() 647 IF (logger%para_env%ionode) THEN 648 unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) 649 ELSE 650 unit_nr = -1 651 ENDIF 652 653 CALL dbcsr_create(matrix, template=matrix_p_in(1), matrix_type=dbcsr_type_no_symmetry) 654 CALL dbcsr_create(matrix_tmp, template=matrix_p_in(1), matrix_type=dbcsr_type_no_symmetry) 655 nspin = SIZE(matrix_p_in) 656 657 DO ispin = 1, nspin 658 step_fac = 1.0_dp 659 frob_norm = 1.0_dp 660 nsave = 0 661 662 CALL dbcsr_copy(matrix_tmp, matrix_p_in(ispin)) 663 CALL dbcsr_copy(matrix_p_out(ispin), matrix_p_in(ispin)) 664! If a BCH history is used make good use of it and do a few steps as a copy and scale update of P 665! else BCH_saved will be 0 and loop is skipped 666 DO i = 1, BCH_saved(ispin) 667 step_fac = step_fac*step_size(ispin) 668 CALL dbcsr_copy(matrix_tmp, matrix_p_out(ispin)) 669 CALL dbcsr_add(matrix_p_out(ispin), matrix_BCH(ispin, i), 1.0_dp, ifac(i)*step_fac) 670 CALL dbcsr_add(matrix_tmp, matrix_p_out(ispin), 1.0_dp, -1.0_dp) 671 frob_norm = dbcsr_frobenius_norm(matrix_tmp) 672 IF (unit_nr .GT. 0) WRITE (unit_nr, "(t3,a,i3,a,f16.8)") "BCH: step", i, " Norm of P_old-Pnew:", frob_norm 673 IF (frob_norm .LT. threshold) EXIT 674 END DO 675 IF (frob_norm .LT. threshold) CYCLE 676 677! If the copy and scale isn't enough compute a few more BCH steps. 20 seems high but except of the first step it will never be close 678 save_BCH = BCH_saved(ispin) == 0 .AND. n_bch_hist .GT. 0 679 DO i = BCH_saved(ispin) + 1, 20 680 step_fac = step_fac*step_size(ispin) 681 !allow for a bit of matrix magic here by exploiting matrix and matrix_tmp 682 !matrix_tmp is alway the previous order of the BCH series 683 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp, matrix_dp(ispin), & 684 0.0_dp, matrix, filter_eps=threshold) 685 686 !(anti)symmetry allows to sum the transposed instead of the full commutator, matrix becomes the latest result 687 688 CALL dbcsr_transposed(matrix_tmp, matrix) 689 CALL dbcsr_add(matrix, matrix_tmp, 1.0_dp, 1.0_dp) 690 691 !Finally, add the new BCH order to P, but store the previous one for a convergence check 692 CALL dbcsr_copy(matrix_tmp, matrix_p_out(ispin)) 693 CALL dbcsr_add(matrix_p_out(ispin), matrix, 1.0_dp, ifac(i)*step_fac) 694 IF (save_BCH .AND. i .LE. n_bch_hist) THEN 695 CALL dbcsr_copy(matrix_BCH(ispin, i), matrix) 696 nsave = i 697 END IF 698 699 CALL dbcsr_add(matrix_tmp, matrix_p_out(ispin), 1.0_dp, -1.0_dp) 700 701 !Stop the BCH-series if two successive P's differ by less the threshold 702 frob_norm = dbcsr_frobenius_norm(matrix_tmp) 703 IF (unit_nr .GT. 0) WRITE (unit_nr, "(t3,a,i3,a,f16.8)") "BCH: step", i, " Norm of P_old-Pnew:", frob_norm 704 IF (frob_norm .LT. threshold) EXIT 705 706 !Copy the latest BCH-matrix on matrix tmp, so we can cycle with all matrices in place 707 CALL dbcsr_copy(matrix_tmp, matrix) 708 CALL dbcsr_filter(matrix_tmp, threshold) 709 END DO 710 BCH_saved(ispin) = nsave 711 IF (unit_nr .GT. 0) WRITE (unit_nr, "(A)") " " 712 END DO 713 714 CALL purify_mcweeny(matrix_p_out, threshold, 1) 715 IF (unit_nr .GT. 0) CALL m_flush(unit_nr) 716 CALL dbcsr_release(matrix_tmp) 717 CALL dbcsr_release(matrix) 718 CALL timestop(handle) 719 END SUBROUTINE update_p_exp 720 721! ************************************************************************************************** 722!> \brief performs a tranformation of a matrix back to/into orthonormal basis 723!> in case of P a scaling of 0.5 has to be applied for closed shell case 724!> \param matrix matrix to be transformed 725!> \param matrix_trafo transformation matrix 726!> \param eps_filter filtering threshold for sparse matrices 727!> \par History 728!> 2012.05 created [Florian Schiffmann] 729!> \author Florian Schiffmann 730! ************************************************************************************************** 731 732 SUBROUTINE transform_matrix_orth(matrix, matrix_trafo, eps_filter) 733 TYPE(dbcsr_type), DIMENSION(:) :: matrix 734 TYPE(dbcsr_type) :: matrix_trafo 735 REAL(KIND=dp) :: eps_filter 736 737 CHARACTER(LEN=*), PARAMETER :: routineN = 'transform_matrix_orth', & 738 routineP = moduleN//':'//routineN 739 740 INTEGER :: handle, ispin 741 TYPE(dbcsr_type) :: matrix_tmp, matrix_work 742 743 CALL timeset(routineN, handle) 744 745 CALL dbcsr_create(matrix_work, template=matrix(1), matrix_type=dbcsr_type_no_symmetry) 746 CALL dbcsr_create(matrix_tmp, template=matrix(1), matrix_type=dbcsr_type_no_symmetry) 747 748 DO ispin = 1, SIZE(matrix) 749 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix(ispin), matrix_trafo, & 750 0.0_dp, matrix_work, filter_eps=eps_filter) 751 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_trafo, matrix_work, & 752 0.0_dp, matrix_tmp, filter_eps=eps_filter) 753 ! symmetrize results (this is again needed to make sure everything is stable) 754 CALL dbcsr_transposed(matrix_work, matrix_tmp) 755 CALL dbcsr_add(matrix_tmp, matrix_work, 0.5_dp, 0.5_dp) 756 CALL dbcsr_copy(matrix(ispin), matrix_tmp) 757 END DO 758 759 CALL dbcsr_release(matrix_tmp) 760 CALL dbcsr_release(matrix_work) 761 CALL timestop(handle) 762 763 END SUBROUTINE 764 765! ************************************************************************************************** 766!> \brief ... 767!> \param curvy_data ... 768! ************************************************************************************************** 769 SUBROUTINE deallocate_curvy_data(curvy_data) 770 TYPE(ls_scf_curvy_type) :: curvy_data 771 772 CHARACTER(LEN=*), PARAMETER :: routineN = 'deallocate_curvy_data', & 773 routineP = moduleN//':'//routineN 774 775 INTEGER :: i, j 776 777 CALL release_dbcsr_array(curvy_data%matrix_dp) 778 CALL release_dbcsr_array(curvy_data%matrix_p) 779 780 IF (ALLOCATED(curvy_data%matrix_psave)) THEN 781 DO i = 1, SIZE(curvy_data%matrix_psave, 1) 782 DO j = 1, 3 783 CALL dbcsr_release(curvy_data%matrix_psave(i, j)) 784 END DO 785 END DO 786 DEALLOCATE (curvy_data%matrix_psave) 787 END IF 788 IF (ALLOCATED(curvy_data%matrix_BCH)) THEN 789 DO i = 1, SIZE(curvy_data%matrix_BCH, 1) 790 DO j = 1, 7 791 CALL dbcsr_release(curvy_data%matrix_BCH(i, j)) 792 END DO 793 END DO 794 DEALLOCATE (curvy_data%matrix_BCH) 795 END IF 796 END SUBROUTINE deallocate_curvy_data 797 798! ************************************************************************************************** 799!> \brief ... 800!> \param matrix ... 801! ************************************************************************************************** 802 SUBROUTINE release_dbcsr_array(matrix) 803 TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: matrix 804 805 CHARACTER(LEN=*), PARAMETER :: routineN = 'release_dbcsr_array', & 806 routineP = moduleN//':'//routineN 807 808 INTEGER :: i 809 810 IF (ALLOCATED(matrix)) THEN 811 DO i = 1, SIZE(matrix) 812 CALL dbcsr_release(matrix(i)) 813 END DO 814 DEALLOCATE (matrix) 815 END IF 816 END SUBROUTINE release_dbcsr_array 817 818! ************************************************************************************************** 819!> \brief ... 820!> \param curvy_data ... 821!> \param matrix_s ... 822!> \param nspins ... 823! ************************************************************************************************** 824 SUBROUTINE init_curvy(curvy_data, matrix_s, nspins) 825 TYPE(ls_scf_curvy_type) :: curvy_data 826 TYPE(dbcsr_type) :: matrix_s 827 INTEGER :: nspins 828 829 CHARACTER(LEN=*), PARAMETER :: routineN = 'init_curvy', routineP = moduleN//':'//routineN 830 831 INTEGER :: ispin, j 832 833 ALLOCATE (curvy_data%matrix_dp(nspins)) 834 ALLOCATE (curvy_data%matrix_p(nspins)) 835 DO ispin = 1, nspins 836 CALL dbcsr_create(curvy_data%matrix_dp(ispin), template=matrix_s, & 837 matrix_type=dbcsr_type_no_symmetry) 838 CALL dbcsr_set(curvy_data%matrix_dp(ispin), 0.0_dp) 839 CALL dbcsr_create(curvy_data%matrix_p(ispin), template=matrix_s, & 840 matrix_type=dbcsr_type_no_symmetry) 841 curvy_data%fix_shift = .FALSE. 842 curvy_data%double_step_size = .TRUE. 843 curvy_data%shift = 1.0_dp 844 curvy_data%BCH_saved = 0 845 curvy_data%step_size = 0.60_dp 846 curvy_data%cg_numer = 0.00_dp 847 curvy_data%cg_denom = 0.00_dp 848 END DO 849 IF (curvy_data%line_search_type == ls_scf_line_search_3point_2d) THEN 850 ALLOCATE (curvy_data%matrix_psave(nspins, 3)) 851 DO ispin = 1, nspins 852 DO j = 1, 3 853 CALL dbcsr_create(curvy_data%matrix_psave(ispin, j), template=matrix_s, & 854 matrix_type=dbcsr_type_no_symmetry) 855 END DO 856 END DO 857 END IF 858 IF (curvy_data%n_bch_hist .GT. 0) THEN 859 ALLOCATE (curvy_data%matrix_BCH(nspins, curvy_data%n_bch_hist)) 860 DO ispin = 1, nspins 861 DO j = 1, curvy_data%n_bch_hist 862 CALL dbcsr_create(curvy_data%matrix_BCH(ispin, j), template=matrix_s, & 863 matrix_type=dbcsr_type_no_symmetry) 864 END DO 865 END DO 866 END IF 867 868 END SUBROUTINE init_curvy 869 870END MODULE dm_ls_scf_curvy 871