1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Routines for Geometry optimization using BFGS algorithm 8! ************************************************************************************************** 9MODULE bfgs_optimizer 10 11 USE atomic_kind_list_types, ONLY: atomic_kind_list_type 12 USE atomic_kind_types, ONLY: get_atomic_kind,& 13 get_atomic_kind_set 14 USE cell_types, ONLY: cell_type,& 15 pbc 16 USE constraint_fxd, ONLY: fix_atom_control 17 USE cp_blacs_env, ONLY: cp_blacs_env_create,& 18 cp_blacs_env_release,& 19 cp_blacs_env_type 20 USE cp_external_control, ONLY: external_control 21 USE cp_files, ONLY: close_file,& 22 open_file 23 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,& 24 cp_fm_transpose 25 USE cp_fm_diag, ONLY: choose_eigv_solver 26 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 27 cp_fm_struct_release,& 28 cp_fm_struct_type 29 USE cp_fm_types, ONLY: cp_fm_create,& 30 cp_fm_get_info,& 31 cp_fm_read_unformatted,& 32 cp_fm_release,& 33 cp_fm_set_all,& 34 cp_fm_to_fm,& 35 cp_fm_type,& 36 cp_fm_write_unformatted 37 USE cp_gemm_interface, ONLY: cp_gemm 38 USE cp_log_handling, ONLY: cp_get_default_logger,& 39 cp_logger_type,& 40 cp_to_string 41 USE cp_output_handling, ONLY: cp_iterate,& 42 cp_p_file,& 43 cp_print_key_finished_output,& 44 cp_print_key_should_output,& 45 cp_print_key_unit_nr 46 USE cp_para_types, ONLY: cp_para_env_type 47 USE cp_subsys_types, ONLY: cp_subsys_get,& 48 cp_subsys_type 49 USE force_env_types, ONLY: force_env_get,& 50 force_env_type 51 USE global_types, ONLY: global_environment_type 52 USE gopt_f_methods, ONLY: gopt_f_ii,& 53 gopt_f_io,& 54 gopt_f_io_finalize,& 55 gopt_f_io_init,& 56 print_geo_opt_header,& 57 print_geo_opt_nc 58 USE gopt_f_types, ONLY: gopt_f_type 59 USE gopt_param_types, ONLY: gopt_param_type 60 USE input_constants, ONLY: default_cell_method_id,& 61 default_ts_method_id 62 USE input_section_types, ONLY: section_vals_get_subs_vals,& 63 section_vals_type,& 64 section_vals_val_get,& 65 section_vals_val_set 66 USE kinds, ONLY: default_path_length,& 67 dp 68 USE machine, ONLY: m_flush,& 69 m_walltime 70 USE message_passing, ONLY: mp_min,& 71 mp_sum 72 USE particle_list_types, ONLY: particle_list_type 73#include "../base/base_uses.f90" 74 75 IMPLICIT NONE 76 PRIVATE 77#include "gopt_f77_methods.h" 78 79 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bfgs_optimizer' 80 LOGICAL, PARAMETER :: debug_this_module = .TRUE. 81 82 PUBLIC :: geoopt_bfgs 83 84CONTAINS 85 86! ************************************************************************************************** 87!> \brief Main driver for BFGS geometry optimizations 88!> \param force_env ... 89!> \param gopt_param ... 90!> \param globenv ... 91!> \param geo_section ... 92!> \param gopt_env ... 93!> \param x0 ... 94! ************************************************************************************************** 95 RECURSIVE SUBROUTINE geoopt_bfgs(force_env, gopt_param, globenv, geo_section, gopt_env, x0) 96 97 TYPE(force_env_type), POINTER :: force_env 98 TYPE(gopt_param_type), POINTER :: gopt_param 99 TYPE(global_environment_type), POINTER :: globenv 100 TYPE(section_vals_type), POINTER :: geo_section 101 TYPE(gopt_f_type), POINTER :: gopt_env 102 REAL(KIND=dp), DIMENSION(:), POINTER :: x0 103 104 CHARACTER(len=*), PARAMETER :: routineN = 'geoopt_bfgs', routineP = moduleN//':'//routineN 105 REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp 106 107 CHARACTER(LEN=5) :: wildcard 108 CHARACTER(LEN=default_path_length) :: hes_filename 109 INTEGER :: handle, hesunit_read, indf, info, & 110 iter_nr, its, maxiter, ndf, nfree, & 111 output_unit 112 LOGICAL :: conv, hesrest, ionode, shell_present, & 113 should_stop, use_mod_hes, use_rfo 114 REAL(KIND=dp) :: ediff, emin, eold, etot, pred, rad, rat, & 115 step, t_diff, t_now, t_old 116 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dg, dr, dx, eigval, gold, work, xold 117 REAL(KIND=dp), DIMENSION(:), POINTER :: g 118 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 119 TYPE(cp_blacs_env_type), POINTER :: blacs_env 120 TYPE(cp_fm_struct_type), POINTER :: fm_struct_hes 121 TYPE(cp_fm_type), POINTER :: eigvec_mat, hess_mat, hess_tmp 122 TYPE(cp_logger_type), POINTER :: logger 123 TYPE(cp_para_env_type), POINTER :: para_env 124 TYPE(cp_subsys_type), POINTER :: subsys 125 TYPE(section_vals_type), POINTER :: print_key, root_section 126 127 NULLIFY (logger, g, blacs_env) 128 logger => cp_get_default_logger() 129 para_env => force_env%para_env 130 root_section => force_env%root_section 131 t_old = m_walltime() 132 133 CALL timeset(routineN, handle) 134 CALL section_vals_val_get(geo_section, "BFGS%TRUST_RADIUS", r_val=rad) 135 print_key => section_vals_get_subs_vals(geo_section, "BFGS%RESTART") 136 ionode = para_env%mepos == para_env%source 137 maxiter = gopt_param%max_iter 138 conv = .FALSE. 139 rat = 0.0_dp 140 wildcard = " BFGS" 141 142 ! Stop if not yet implemented 143 SELECT CASE (gopt_env%type_id) 144 CASE (default_ts_method_id) 145 CPABORT("BFGS method not yet working with DIMER") 146 END SELECT 147 148 CALL section_vals_val_get(geo_section, "BFGS%USE_RAT_FUN_OPT", l_val=use_rfo) 149 CALL section_vals_val_get(geo_section, "BFGS%USE_MODEL_HESSIAN", l_val=use_mod_hes) 150 CALL section_vals_val_get(geo_section, "BFGS%RESTART_HESSIAN", l_val=hesrest) 151 output_unit = cp_print_key_unit_nr(logger, geo_section, "PRINT%PROGRAM_RUN_INFO", & 152 extension=".geoLog") 153 IF (output_unit > 0) THEN 154 IF (use_rfo) THEN 155 WRITE (UNIT=output_unit, FMT="(/,T2,A,T78,A3)") & 156 "BFGS| Use rational function optimization for step estimation: ", "YES" 157 ELSE 158 WRITE (UNIT=output_unit, FMT="(/,T2,A,T78,A3)") & 159 "BFGS| Use rational function optimization for step estimation: ", " NO" 160 END IF 161 IF (use_mod_hes) THEN 162 WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") & 163 "BFGS| Use model Hessian for initial guess: ", "YES" 164 ELSE 165 WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") & 166 "BFGS| Use model Hessian for initial guess: ", " NO" 167 END IF 168 IF (hesrest) THEN 169 WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") & 170 "BFGS| Restart Hessian: ", "YES" 171 ELSE 172 WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") & 173 "BFGS| Restart Hessian: ", " NO" 174 END IF 175 WRITE (UNIT=output_unit, FMT="(T2,A,T61,F20.3)") & 176 "BFGS| Trust radius: ", rad 177 END IF 178 179 ndf = SIZE(x0) 180 nfree = gopt_env%nfree 181 IF (ndf > 3000) & 182 CALL cp_warn(__LOCATION__, & 183 "The dimension of the Hessian matrix ("// & 184 TRIM(ADJUSTL(cp_to_string(ndf)))//") is greater than 3000. "// & 185 "The diagonalisation of the full Hessian matrix needed for BFGS "// & 186 "is computationally expensive. You should consider to use the linear "// & 187 "scaling variant L-BFGS instead.") 188 189 ! Initialize hessian (hes = unitary matrix or model hessian ) 190 CALL cp_blacs_env_create(blacs_env, para_env, globenv%blacs_grid_layout, & 191 globenv%blacs_repeatable) 192 CALL cp_fm_struct_create(fm_struct_hes, para_env=para_env, context=blacs_env, & 193 nrow_global=ndf, ncol_global=ndf) 194 CALL cp_fm_create(hess_mat, fm_struct_hes, name="hess_mat") 195 CALL cp_fm_create(hess_tmp, fm_struct_hes, name="hess_tmp") 196 CALL cp_fm_create(eigvec_mat, fm_struct_hes, name="eigvec_mat") 197 ALLOCATE (eigval(ndf)) 198 eigval(:) = 0.0_dp 199 200 CALL force_env_get(force_env=force_env, subsys=subsys) 201 CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds) 202 CALL get_atomic_kind_set(atomic_kind_set=atomic_kinds%els, shell_present=shell_present) 203 IF (use_mod_hes) THEN 204 IF (shell_present) THEN 205 CALL cp_warn(__LOCATION__, & 206 "No model Hessian is available for core-shell models. "// & 207 "A unit matrix is used as the initial Hessian.") 208 use_mod_hes = .FALSE. 209 END IF 210 IF (gopt_env%type_id == default_cell_method_id) THEN 211 CALL cp_warn(__LOCATION__, & 212 "No model Hessian is available for cell optimizations. "// & 213 "A unit matrix is used as the initial Hessian.") 214 use_mod_hes = .FALSE. 215 END IF 216 END IF 217 218 IF (use_mod_hes) THEN 219 CALL cp_fm_set_all(hess_mat, alpha=zero, beta=0.00_dp) 220 CALL construct_initial_hess(gopt_env%force_env, hess_mat) 221 CALL cp_fm_to_fm(hess_mat, hess_tmp) 222 CALL choose_eigv_solver(hess_tmp, eigvec_mat, eigval, info=info) 223 ! In rare cases the diagonalization of hess_mat fails (bug in scalapack?) 224 IF (info /= 0) THEN 225 CALL cp_fm_set_all(hess_mat, alpha=zero, beta=one) 226 IF (output_unit > 0) WRITE (output_unit, *) & 227 "BFGS: Matrix diagonalization failed, using unity as model Hessian." 228 ELSE 229 DO its = 1, SIZE(eigval) 230 IF (eigval(its) .LT. 0.1_dp) eigval(its) = 0.1_dp 231 END DO 232 CALL cp_fm_to_fm(eigvec_mat, hess_tmp) 233 CALL cp_fm_column_scale(eigvec_mat, eigval) 234 CALL cp_gemm("N", "T", ndf, ndf, ndf, one, hess_tmp, eigvec_mat, zero, hess_mat) 235 ENDIF 236 ELSE 237 CALL cp_fm_set_all(hess_mat, alpha=zero, beta=one) 238 END IF 239 240 ALLOCATE (xold(ndf)) 241 xold(:) = x0(:) 242 243 ALLOCATE (g(ndf)) 244 g(:) = 0.0_dp 245 246 ALLOCATE (gold(ndf)) 247 gold(:) = 0.0_dp 248 249 ALLOCATE (dx(ndf)) 250 dx(:) = 0.0_dp 251 252 ALLOCATE (dg(ndf)) 253 dg(:) = 0.0_dp 254 255 ALLOCATE (work(ndf)) 256 work(:) = 0.0_dp 257 258 ALLOCATE (dr(ndf)) 259 dr(:) = 0.0_dp 260 261 ! Geometry optimization starts now 262 CALL cp_iterate(logger%iter_info, increment=0, iter_nr_out=iter_nr) 263 CALL print_geo_opt_header(gopt_env, output_unit, wildcard) 264 265 ! Calculate Energy & Gradients 266 CALL cp_eval_at(gopt_env, x0, etot, g, gopt_env%force_env%para_env%mepos, & 267 .FALSE., gopt_env%force_env%para_env) 268 269 ! Print info at time 0 270 emin = etot 271 t_now = m_walltime() 272 t_diff = t_now - t_old 273 t_old = t_now 274 CALL gopt_f_io_init(gopt_env, output_unit, etot, wildcard=wildcard, its=iter_nr, used_time=t_diff) 275 DO its = iter_nr + 1, maxiter 276 CALL cp_iterate(logger%iter_info, last=(its == maxiter)) 277 CALL section_vals_val_set(geo_section, "STEP_START_VAL", i_val=its) 278 CALL gopt_f_ii(its, output_unit) 279 280 ! Hessian update/restarting 281 IF (((its - iter_nr) == 1) .AND. hesrest) THEN 282 IF (ionode) THEN 283 CALL section_vals_val_get(geo_section, "BFGS%RESTART_FILE_NAME", c_val=hes_filename) 284 CALL open_file(file_name=hes_filename, file_status="OLD", & 285 file_form="UNFORMATTED", file_action="READ", unit_number=hesunit_read) 286 END IF 287 CALL cp_fm_read_unformatted(hess_mat, hesunit_read) 288 IF (ionode) CALL close_file(unit_number=hesunit_read) 289 ELSE 290 IF ((its - iter_nr) > 1) THEN 291 DO indf = 1, ndf 292 dx(indf) = x0(indf) - xold(indf) 293 dg(indf) = g(indf) - gold(indf) 294 END DO 295 296 CALL bfgs(ndf, dx, dg, hess_mat, work, para_env) 297 298 !Possibly dump the Hessian file 299 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 300 CALL write_bfgs_hessian(geo_section, hess_mat, logger) 301 ENDIF 302 ENDIF 303 END IF 304 305 ! Setting the present positions & gradients as old 306 xold(:) = x0 307 gold(:) = g 308 309 ! Copying hessian hes to (ndf x ndf) matrix hes_mat for diagonalization 310 CALL cp_fm_to_fm(hess_mat, hess_tmp) 311 312 CALL choose_eigv_solver(hess_tmp, eigvec_mat, eigval, info=info) 313 314 ! In rare cases the diagonalization of hess_mat fails (bug in scalapack?) 315 IF (info /= 0) THEN 316 IF (output_unit > 0) WRITE (output_unit, *) & 317 "BFGS: Matrix diagonalization failed, resetting Hessian to unity." 318 CALL cp_fm_set_all(hess_mat, alpha=zero, beta=one) 319 CALL cp_fm_to_fm(hess_mat, hess_tmp) 320 CALL choose_eigv_solver(hess_tmp, eigvec_mat, eigval) 321 END IF 322 323 IF (use_rfo) THEN 324 CALL set_hes_eig(ndf, eigval, work) 325 dx(:) = eigval 326 CALL rat_fun_opt(ndf, dg, eigval, work, eigvec_mat, g, para_env) 327 END IF 328 CALL geoopt_get_step(ndf, eigval, eigvec_mat, hess_tmp, dr, g, para_env, use_rfo) 329 CALL trust_radius(ndf, step, rad, rat, dr, output_unit) 330 331 ! Update the atomic positions 332 x0 = x0 + dr 333 334 CALL energy_predict(ndf, work, hess_mat, dr, g, conv, pred, para_env) 335 eold = etot 336 337 ! Energy & Gradients at new step 338 CALL cp_eval_at(gopt_env, x0, etot, g, gopt_env%force_env%para_env%mepos, & 339 .FALSE., gopt_env%force_env%para_env) 340 ediff = etot - eold 341 342 ! check for an external exit command 343 CALL external_control(should_stop, "GEO", globenv=globenv) 344 IF (should_stop) EXIT 345 346 ! Some IO and Convergence check 347 t_now = m_walltime() 348 t_diff = t_now - t_old 349 t_old = t_now 350 CALL gopt_f_io(gopt_env, force_env, root_section, its, etot, output_unit, & 351 eold, emin, wildcard, gopt_param, ndf, dr, g, conv, pred, rat, & 352 step, rad, used_time=t_diff) 353 354 IF (conv .OR. (its == maxiter)) EXIT 355 IF (etot < emin) emin = etot 356 IF (use_rfo) CALL update_trust_rad(rat, rad, step, ediff) 357 END DO 358 359 IF (its == maxiter .AND. (.NOT. conv)) THEN 360 CALL print_geo_opt_nc(gopt_env, output_unit) 361 END IF 362 363 ! Write final information, if converged 364 CALL cp_iterate(logger%iter_info, last=.TRUE., increment=0) 365 CALL write_bfgs_hessian(geo_section, hess_mat, logger) 366 CALL gopt_f_io_finalize(gopt_env, force_env, x0, conv, its, root_section, & 367 gopt_env%force_env%para_env, gopt_env%force_env%para_env%mepos, output_unit) 368 369 CALL cp_fm_struct_release(fm_struct_hes) 370 CALL cp_fm_release(hess_mat) 371 CALL cp_fm_release(eigvec_mat) 372 CALL cp_fm_release(hess_tmp) 373 374 CALL cp_blacs_env_release(blacs_env) 375 DEALLOCATE (xold) 376 DEALLOCATE (g) 377 DEALLOCATE (gold) 378 DEALLOCATE (dx) 379 DEALLOCATE (dg) 380 DEALLOCATE (eigval) 381 DEALLOCATE (work) 382 DEALLOCATE (dr) 383 384 CALL cp_print_key_finished_output(output_unit, logger, geo_section, & 385 "PRINT%PROGRAM_RUN_INFO") 386 CALL timestop(handle) 387 388 END SUBROUTINE geoopt_bfgs 389 390! ************************************************************************************************** 391!> \brief ... 392!> \param ndf ... 393!> \param dg ... 394!> \param eigval ... 395!> \param work ... 396!> \param eigvec_mat ... 397!> \param g ... 398!> \param para_env ... 399! ************************************************************************************************** 400 SUBROUTINE rat_fun_opt(ndf, dg, eigval, work, eigvec_mat, g, para_env) 401 402 INTEGER, INTENT(IN) :: ndf 403 REAL(KIND=dp), INTENT(INOUT) :: dg(ndf), eigval(ndf), work(ndf) 404 TYPE(cp_fm_type), POINTER :: eigvec_mat 405 REAL(KIND=dp), INTENT(INOUT) :: g(ndf) 406 TYPE(cp_para_env_type), OPTIONAL, POINTER :: para_env 407 408 CHARACTER(LEN=*), PARAMETER :: routineN = 'rat_fun_opt', routineP = moduleN//':'//routineN 409 REAL(KIND=dp), PARAMETER :: one = 1.0_dp 410 411 INTEGER :: handle, i, indf, iref, iter, j, k, l, & 412 maxit, ncol_local, nrow_local 413 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 414 LOGICAL :: bisec, fail, set 415 REAL(KIND=dp) :: fun, fun1, fun2, fun3, fung, lam1, lam2, & 416 ln, lp, ssize, step, stol 417 REAL(KIND=dp), DIMENSION(:, :), POINTER :: local_data 418 419 CALL timeset(routineN, handle) 420 421 stol = 1.0E-8_dp 422 ssize = 0.2_dp 423 maxit = 999 424 fail = .FALSE. 425 bisec = .FALSE. 426 427 dg = 0._dp 428 429 CALL cp_fm_get_info(eigvec_mat, row_indices=row_indices, col_indices=col_indices, & 430 local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local) 431 432 DO i = 1, nrow_local 433 j = row_indices(i) 434 DO k = 1, ncol_local 435 l = col_indices(k) 436 dg(l) = dg(l) + local_data(i, k)*g(j) 437 END DO 438 END DO 439 CALL mp_sum(dg, para_env%group) 440 441 set = .FALSE. 442 443 DO 444 445! calculating Lamda 446 447 lp = 0.0_dp 448 iref = 1 449 ln = 0.0_dp 450 IF (eigval(iref) < 0.0_dp) ln = eigval(iref) - 0.01_dp 451 452 iter = 0 453 DO 454 iter = iter + 1 455 fun = 0.0_dp 456 fung = 0.0_dp 457 DO indf = 1, ndf 458 fun = fun + dg(indf)**2/(ln - eigval(indf)) 459 fung = fung - dg(indf)**2/(ln - eigval(indf)**2) 460 END DO 461 fun = fun - ln 462 fung = fung - one 463 step = fun/fung 464 ln = ln - step 465 IF (ABS(step) < stol) GOTO 200 466 IF (iter >= maxit) EXIT 467 END DO 468100 CONTINUE 469 bisec = .TRUE. 470 iter = 0 471 maxit = 9999 472 lam1 = 0.0_dp 473 IF (eigval(iref) < 0.0_dp) lam1 = eigval(iref) - 0.01_dp 474 fun1 = 0.0_dp 475 DO indf = 1, ndf 476 fun1 = fun1 + dg(indf)**2/(lam1 - eigval(indf)) 477 END DO 478 fun1 = fun1 - lam1 479 step = ABS(lam1)/1000.0_dp 480 IF (step < ssize) step = ssize 481 DO 482 iter = iter + 1 483 IF (iter > maxit) THEN 484 ln = 0.0_dp 485 lp = 0.0_dp 486 fail = .TRUE. 487 GOTO 300 488 END IF 489 fun2 = 0.0_dp 490 lam2 = lam1 - iter*step 491 DO indf = 1, ndf 492 fun2 = fun2 + eigval(indf)**2/(lam2 - eigval(indf)) 493 END DO 494 fun2 = fun2 - lam2 495 IF (fun2*fun1 < 0.0_dp) THEN 496 iter = 0 497 DO 498 iter = iter + 1 499 IF (iter > maxit) THEN 500 ln = 0.0_dp 501 lp = 0.0_dp 502 fail = .TRUE. 503 GOTO 300 504 END IF 505 step = (lam1 + lam2)/2 506 fun3 = 0.0_dp 507 DO indf = 1, ndf 508 fun3 = fun3 + dg(indf)**2/(step - eigval(indf)) 509 END DO 510 fun3 = fun3 - step 511 512 IF (ABS(step - lam2) < stol) THEN 513 ln = step 514 GOTO 200 515 END IF 516 517 IF (fun3*fun1 < stol) THEN 518 lam2 = step 519 ELSE 520 lam1 = step 521 END IF 522 END DO 523 END IF 524 END DO 525 526200 CONTINUE 527 IF ((ln > eigval(iref)) .OR. ((ln > 0.0_dp) .AND. & 528 (eigval(iref) > 0.0_dp))) THEN 529 530 IF (.NOT. bisec) GOTO 100 531 ln = 0.0_dp 532 lp = 0.0_dp 533 fail = .TRUE. 534 END IF 535 536300 CONTINUE 537 538 IF (fail .AND. .NOT. set) THEN 539 set = .TRUE. 540 DO indf = 1, ndf 541 eigval(indf) = eigval(indf)*work(indf) 542 END DO 543 CYCLE 544 END IF 545 546 IF (.NOT. set) THEN 547 work(1:ndf) = one 548 ENDIF 549 550 DO indf = 1, ndf 551 eigval(indf) = eigval(indf) - ln 552 END DO 553 EXIT 554 END DO 555 556 CALL timestop(handle) 557 558 END SUBROUTINE rat_fun_opt 559 560! ************************************************************************************************** 561!> \brief ... 562!> \param ndf ... 563!> \param dx ... 564!> \param dg ... 565!> \param hess_mat ... 566!> \param work ... 567!> \param para_env ... 568! ************************************************************************************************** 569 SUBROUTINE bfgs(ndf, dx, dg, hess_mat, work, para_env) 570 INTEGER, INTENT(IN) :: ndf 571 REAL(KIND=dp), INTENT(INOUT) :: dx(ndf), dg(ndf) 572 TYPE(cp_fm_type), POINTER :: hess_mat 573 REAL(KIND=dp), INTENT(INOUT) :: work(ndf) 574 TYPE(cp_para_env_type), OPTIONAL, POINTER :: para_env 575 576 CHARACTER(LEN=*), PARAMETER :: routineN = 'bfgs', routineP = moduleN//':'//routineN 577 REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp 578 579 INTEGER :: handle, i, j, k, l, ncol_local, & 580 nrow_local 581 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 582 REAL(KIND=dp) :: DDOT, dxw, gdx 583 REAL(KIND=dp), DIMENSION(:, :), POINTER :: local_hes 584 585 CALL timeset(routineN, handle) 586 587 CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, & 588 local_data=local_hes, nrow_local=nrow_local, ncol_local=ncol_local) 589 590 work = zero 591 DO i = 1, nrow_local 592 j = row_indices(i) 593 DO k = 1, ncol_local 594 l = col_indices(k) 595 work(j) = work(j) + local_hes(i, k)*dx(l) 596 END DO 597 END DO 598 599 CALL mp_sum(work, para_env%group) 600 601 gdx = DDOT(ndf, dg, 1, dx, 1) 602 gdx = one/gdx 603 dxw = DDOT(ndf, dx, 1, work, 1) 604 dxw = one/dxw 605 606 DO i = 1, nrow_local 607 j = row_indices(i) 608 DO k = 1, ncol_local 609 l = col_indices(k) 610 local_hes(i, k) = local_hes(i, k) + gdx*dg(j)*dg(l) - & 611 dxw*work(j)*work(l) 612 END DO 613 END DO 614 615 CALL timestop(handle) 616 617 END SUBROUTINE bfgs 618 619! ************************************************************************************************** 620!> \brief ... 621!> \param ndf ... 622!> \param eigval ... 623!> \param work ... 624! ************************************************************************************************** 625 SUBROUTINE set_hes_eig(ndf, eigval, work) 626 INTEGER, INTENT(IN) :: ndf 627 REAL(KIND=dp), INTENT(INOUT) :: eigval(ndf), work(ndf) 628 629 CHARACTER(LEN=*), PARAMETER :: routineN = 'set_hes_eig', routineP = moduleN//':'//routineN 630 REAL(KIND=dp), PARAMETER :: max_neg = -0.5_dp, max_pos = 5.0_dp, & 631 min_eig = 0.005_dp, one = 1.0_dp 632 633 INTEGER :: handle, indf 634 LOGICAL :: neg 635 636 CALL timeset(routineN, handle) 637 638 DO indf = 1, ndf 639 IF (eigval(indf) < 0.0_dp) neg = .TRUE. 640 IF (eigval(indf) > 1000.0_dp) eigval(indf) = 1000.0_dp 641 END DO 642 DO indf = 1, ndf 643 IF (eigval(indf) < 0.0_dp) THEN 644 IF (eigval(indf) < max_neg) THEN 645 eigval(indf) = max_neg 646 ELSE IF (eigval(indf) > -min_eig) THEN 647 eigval(indf) = -min_eig 648 END IF 649 ELSE IF (eigval(indf) < 1000.0_dp) THEN 650 IF (eigval(indf) < min_eig) THEN 651 eigval(indf) = min_eig 652 ELSE IF (eigval(indf) > max_pos) THEN 653 eigval(indf) = max_pos 654 END IF 655 END IF 656 END DO 657 658 DO indf = 1, ndf 659 IF (eigval(indf) < 0.0_dp) THEN 660 work(indf) = -one 661 ELSE 662 work(indf) = one 663 END IF 664 END DO 665 666 CALL timestop(handle) 667 668 END SUBROUTINE set_hes_eig 669 670! ************************************************************************************************** 671!> \brief ... 672!> \param ndf ... 673!> \param eigval ... 674!> \param eigvec_mat ... 675!> \param hess_tmp ... 676!> \param dr ... 677!> \param g ... 678!> \param para_env ... 679!> \param use_rfo ... 680! ************************************************************************************************** 681 SUBROUTINE geoopt_get_step(ndf, eigval, eigvec_mat, hess_tmp, dr, g, para_env, use_rfo) 682 683 INTEGER, INTENT(IN) :: ndf 684 REAL(KIND=dp), INTENT(INOUT) :: eigval(ndf) 685 TYPE(cp_fm_type), POINTER :: eigvec_mat, hess_tmp 686 REAL(KIND=dp), INTENT(INOUT) :: dr(ndf), g(ndf) 687 TYPE(cp_para_env_type), OPTIONAL, POINTER :: para_env 688 LOGICAL :: use_rfo 689 690 REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp 691 692 INTEGER :: i, indf, j, k, l, ncol_local, nrow_local 693 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 694 REAL(KIND=dp), DIMENSION(:, :), POINTER :: local_data 695 TYPE(cp_fm_struct_type), POINTER :: matrix_struct 696 TYPE(cp_fm_type), POINTER :: tmp 697 698 CALL cp_fm_to_fm(eigvec_mat, hess_tmp) 699 IF (use_rfo) THEN 700 DO indf = 1, ndf 701 eigval(indf) = one/eigval(indf) 702 END DO 703 ELSE 704 DO indf = 1, ndf 705 eigval(indf) = one/MAX(0.0001_dp, eigval(indf)) 706 END DO 707 END IF 708 709 CALL cp_fm_column_scale(hess_tmp, eigval) 710 CALL cp_fm_get_info(eigvec_mat, matrix_struct=matrix_struct) 711 CALL cp_fm_create(tmp, matrix_struct, name="tmp") 712 713 CALL cp_gemm("N", "T", ndf, ndf, ndf, one, hess_tmp, eigvec_mat, zero, tmp) 714 715 CALL cp_fm_transpose(tmp, hess_tmp) 716 CALL cp_fm_release(tmp) 717 718! ** New step ** 719 720 CALL cp_fm_get_info(hess_tmp, row_indices=row_indices, col_indices=col_indices, & 721 local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local) 722 723 dr = 0.0_dp 724 DO i = 1, nrow_local 725 j = row_indices(i) 726 DO k = 1, ncol_local 727 l = col_indices(k) 728 dr(j) = dr(j) - local_data(i, k)*g(l) 729 END DO 730 END DO 731 732 CALL mp_sum(dr, para_env%group) 733 734 END SUBROUTINE geoopt_get_step 735 736! ************************************************************************************************** 737!> \brief ... 738!> \param ndf ... 739!> \param step ... 740!> \param rad ... 741!> \param rat ... 742!> \param dr ... 743!> \param output_unit ... 744! ************************************************************************************************** 745 SUBROUTINE trust_radius(ndf, step, rad, rat, dr, output_unit) 746 INTEGER, INTENT(IN) :: ndf 747 REAL(KIND=dp), INTENT(INOUT) :: step, rad, rat, dr(ndf) 748 INTEGER, INTENT(IN) :: output_unit 749 750 CHARACTER(LEN=*), PARAMETER :: routineN = 'trust_radius', routineP = moduleN//':'//routineN 751 REAL(KIND=dp), PARAMETER :: one = 1.0_dp 752 753 INTEGER :: handle 754 REAL(KIND=dp) :: scal 755 756 CALL timeset(routineN, handle) 757 758 step = MAXVAL(ABS(dr)) 759 scal = MAX(one, rad/step) 760 761 IF (step > rad) THEN 762 rat = rad/step 763 CALL DSCAL(ndf, rat, dr, 1) 764 step = rad 765 IF (output_unit > 0) THEN 766 WRITE (unit=output_unit, FMT="(/,T2,A,F8.5)") & 767 " Step is scaled; Scaling factor = ", rat 768 CALL m_flush(output_unit) 769 ENDIF 770 END IF 771 CALL timestop(handle) 772 773 END SUBROUTINE trust_radius 774 775! ************************************************************************************************** 776!> \brief ... 777!> \param ndf ... 778!> \param work ... 779!> \param hess_mat ... 780!> \param dr ... 781!> \param g ... 782!> \param conv ... 783!> \param pred ... 784!> \param para_env ... 785! ************************************************************************************************** 786 SUBROUTINE energy_predict(ndf, work, hess_mat, dr, g, conv, pred, para_env) 787 788 INTEGER, INTENT(IN) :: ndf 789 REAL(KIND=dp), INTENT(INOUT) :: work(ndf) 790 TYPE(cp_fm_type), POINTER :: hess_mat 791 REAL(KIND=dp), INTENT(INOUT) :: dr(ndf), g(ndf) 792 LOGICAL, INTENT(INOUT) :: conv 793 REAL(KIND=dp), INTENT(INOUT) :: pred 794 TYPE(cp_para_env_type), POINTER :: para_env 795 796 CHARACTER(LEN=*), PARAMETER :: routineN = 'energy_predict', routineP = moduleN//':'//routineN 797 REAL(KIND=dp), PARAMETER :: zero = 0.0_dp 798 799 INTEGER :: handle, i, j, k, l, ncol_local, & 800 nrow_local 801 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 802 REAL(KIND=dp) :: DDOT, ener1, ener2 803 REAL(KIND=dp), DIMENSION(:, :), POINTER :: local_data 804 805 CALL timeset(routineN, handle) 806 807 ener1 = DDOT(ndf, g, 1, dr, 1) 808 809 CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, & 810 local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local) 811 812 work = zero 813 DO i = 1, nrow_local 814 j = row_indices(i) 815 DO k = 1, ncol_local 816 l = col_indices(k) 817 work(j) = work(j) + local_data(i, k)*dr(l) 818 END DO 819 END DO 820 821 CALL mp_sum(work, para_env%group) 822 ener2 = DDOT(ndf, dr, 1, work, 1) 823 pred = ener1 + 0.5_dp*ener2 824 conv = .FALSE. 825 CALL timestop(handle) 826 827 END SUBROUTINE energy_predict 828 829! ************************************************************************************************** 830!> \brief ... 831!> \param rat ... 832!> \param rad ... 833!> \param step ... 834!> \param ediff ... 835! ************************************************************************************************** 836 SUBROUTINE update_trust_rad(rat, rad, step, ediff) 837 838 REAL(KIND=dp), INTENT(INOUT) :: rat, rad, step, ediff 839 840 CHARACTER(LEN=*), PARAMETER :: routineN = 'update_trust_rad', & 841 routineP = moduleN//':'//routineN 842 REAL(KIND=dp), PARAMETER :: max_trust = 1.0_dp, min_trust = 0.1_dp 843 844 INTEGER :: handle 845 846 CALL timeset(routineN, handle) 847 848 IF (rat > 4.0_dp) THEN 849 IF (ediff < 0.0_dp) THEN 850 rad = step*0.5_dp 851 ELSE 852 rad = step*0.25_dp 853 END IF 854 ELSE IF (rat > 2.0_dp) THEN 855 IF (ediff < 0.0_dp) THEN 856 rad = step*0.75_dp 857 ELSE 858 rad = step*0.5_dp 859 END IF 860 ELSE IF (rat > 4.0_dp/3.0_dp) THEN 861 IF (ediff < 0.0_dp) THEN 862 rad = step 863 ELSE 864 rad = step*0.75_dp 865 END IF 866 ELSE IF (rat > 10.0_dp/9.0_dp) THEN 867 IF (ediff < 0.0_dp) THEN 868 rad = step*1.25_dp 869 ELSE 870 rad = step 871 END IF 872 ELSE IF (rat > 0.9_dp) THEN 873 IF (ediff < 0.0_dp) THEN 874 rad = step*1.5_dp 875 ELSE 876 rad = step*1.25_dp 877 END IF 878 ELSE IF (rat > 0.75_dp) THEN 879 IF (ediff < 0.0_dp) THEN 880 rad = step*1.25_dp 881 ELSE 882 rad = step 883 END IF 884 ELSE IF (rat > 0.5_dp) THEN 885 IF (ediff < 0.0_dp) THEN 886 rad = step 887 ELSE 888 rad = step*0.75_dp 889 END IF 890 ELSE IF (rat > 0.25_dp) THEN 891 IF (ediff < 0.0_dp) THEN 892 rad = step*0.75_dp 893 ELSE 894 rad = step*0.5_dp 895 END IF 896 ELSE IF (ediff < 0.0_dp) THEN 897 rad = step*0.5_dp 898 ELSE 899 rad = step*0.25_dp 900 END IF 901 902 rad = MAX(rad, min_trust) 903 rad = MIN(rad, max_trust) 904 CALL timestop(handle) 905 906 END SUBROUTINE update_trust_rad 907 908! ************************************************************************************************** 909 910! ************************************************************************************************** 911!> \brief ... 912!> \param geo_section ... 913!> \param hess_mat ... 914!> \param logger ... 915! ************************************************************************************************** 916 SUBROUTINE write_bfgs_hessian(geo_section, hess_mat, logger) 917 918 TYPE(section_vals_type), POINTER :: geo_section 919 TYPE(cp_fm_type), POINTER :: hess_mat 920 TYPE(cp_logger_type), POINTER :: logger 921 922 CHARACTER(LEN=*), PARAMETER :: routineN = 'write_bfgs_hessian', & 923 routineP = moduleN//':'//routineN 924 925 INTEGER :: handle, hesunit 926 927 CALL timeset(routineN, handle) 928 929 hesunit = cp_print_key_unit_nr(logger, geo_section, "BFGS%RESTART", & 930 extension=".Hessian", file_form="UNFORMATTED", file_action="WRITE", & 931 file_position="REWIND") 932 933 CALL cp_fm_write_unformatted(hess_mat, hesunit) 934 935 CALL cp_print_key_finished_output(hesunit, logger, geo_section, "BFGS%RESTART") 936 937 CALL timestop(handle) 938 939 END SUBROUTINE write_bfgs_hessian 940! ************************************************************************************************** 941 942! ************************************************************************************************** 943!> \brief ... 944!> \param force_env ... 945!> \param hess_mat ... 946! ************************************************************************************************** 947 SUBROUTINE construct_initial_hess(force_env, hess_mat) 948 949 TYPE(force_env_type), POINTER :: force_env 950 TYPE(cp_fm_type), POINTER :: hess_mat 951 952 CHARACTER(LEN=*), PARAMETER :: routineN = 'construct_initial_hess', & 953 routineP = moduleN//':'//routineN 954 955 INTEGER :: i, iat_col, iat_row, iglobal, iind, j, & 956 jat_row, jglobal, jind, k, natom, & 957 ncol_local, nrow_local, z 958 INTEGER, ALLOCATABLE, DIMENSION(:) :: at_row 959 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 960 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: d_ij, rho_ij 961 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: r_ij 962 REAL(KIND=dp), DIMENSION(3, 3) :: alpha, r0 963 REAL(KIND=dp), DIMENSION(:, :), POINTER :: fixed, local_data 964 TYPE(cell_type), POINTER :: cell 965 TYPE(cp_subsys_type), POINTER :: subsys 966 TYPE(particle_list_type), POINTER :: particles 967 968 CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell) 969 CALL cp_subsys_get(subsys, & 970 particles=particles) 971 972 alpha(1, :) = (/1._dp, 0.3949_dp, 0.3949_dp/) 973 alpha(2, :) = (/0.3494_dp, 0.2800_dp, 0.2800_dp/) 974 alpha(3, :) = (/0.3494_dp, 0.2800_dp, 0.1800_dp/) 975 976 r0(1, :) = (/1.35_dp, 2.10_dp, 2.53_dp/) 977 r0(2, :) = (/2.10_dp, 2.87_dp, 3.40_dp/) 978 r0(3, :) = (/2.53_dp, 3.40_dp, 3.40_dp/) 979 980 CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, & 981 local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local) 982 natom = particles%n_els 983 ALLOCATE (at_row(natom)) 984 ALLOCATE (rho_ij(natom, natom)) 985 ALLOCATE (d_ij(natom, natom)) 986 ALLOCATE (r_ij(natom, natom, 3)) 987 ALLOCATE (fixed(3, natom)) 988 fixed = 1.0_dp 989 CALL fix_atom_control(force_env, fixed) 990 DO i = 1, 3 991 CALL mp_min(fixed(i, :), hess_mat%matrix_struct%para_env%group) 992 END DO 993 rho_ij = 0 994 !XXXX insert proper rows !XXX 995 at_row = 3 996 DO i = 1, natom 997 CALL get_atomic_kind(atomic_kind=particles%els(i)%atomic_kind, z=z) 998 IF (z .LE. 10) at_row(i) = 2 999 IF (z .LE. 2) at_row(i) = 1 1000 END DO 1001 DO i = 2, natom 1002 iat_row = at_row(i) 1003 DO j = 1, i - 1 1004 jat_row = at_row(j) 1005 !pbc for a distance vector 1006 r_ij(j, i, :) = pbc(particles%els(i)%r, particles%els(j)%r, cell) 1007 r_ij(i, j, :) = -r_ij(j, i, :) 1008 d_ij(j, i) = SQRT(DOT_PRODUCT(r_ij(j, i, :), r_ij(j, i, :))) 1009 d_ij(i, j) = d_ij(j, i) 1010 rho_ij(j, i) = EXP(alpha(jat_row, iat_row)*(r0(jat_row, iat_row)**2 - d_ij(j, i)**2)) 1011 rho_ij(i, j) = rho_ij(j, i) 1012 END DO 1013 END DO 1014 DO i = 1, ncol_local 1015 iglobal = col_indices(i) 1016 iind = MOD(iglobal - 1, 3) + 1 1017 iat_col = (iglobal + 2)/3 1018 IF (iat_col .GT. natom) CYCLE 1019 DO j = 1, nrow_local 1020 jglobal = row_indices(j) 1021 jind = MOD(jglobal - 1, 3) + 1 1022 iat_row = (jglobal + 2)/3 1023 IF (iat_row .GT. natom) CYCLE 1024 IF (iat_row .NE. iat_col) THEN 1025 IF (d_ij(iat_row, iat_col) .LT. 6.0_dp) & 1026 local_data(j, i) = local_data(j, i) + & 1027 angle_second_deriv(r_ij, d_ij, rho_ij, iind, jind, iat_col, iat_row, natom) 1028 ELSE 1029 local_data(j, i) = local_data(j, i) + & 1030 angle_second_deriv(r_ij, d_ij, rho_ij, iind, jind, iat_col, iat_row, natom) 1031 END IF 1032 IF (iat_col .NE. iat_row) THEN 1033 IF (d_ij(iat_row, iat_col) .LT. 6.0_dp) & 1034 local_data(j, i) = local_data(j, i) - & 1035 dist_second_deriv(r_ij(iat_col, iat_row, :), & 1036 iind, jind, d_ij(iat_row, iat_col), rho_ij(iat_row, iat_col)) 1037 ELSE 1038 DO k = 1, natom 1039 IF (k == iat_col) CYCLE 1040 IF (d_ij(iat_row, k) .LT. 6.0_dp) & 1041 local_data(j, i) = local_data(j, i) + & 1042 dist_second_deriv(r_ij(iat_col, k, :), & 1043 iind, jind, d_ij(iat_row, k), rho_ij(iat_row, k)) 1044 END DO 1045 END IF 1046 IF (fixed(jind, iat_row) .LT. 0.5_dp .OR. fixed(iind, iat_col) .LT. 0.5_dp) THEN 1047 local_data(j, i) = 0.0_dp 1048 IF (jind == iind .AND. iat_row == iat_col) local_data(j, i) = 1.0_dp 1049 END IF 1050 END DO 1051 END DO 1052 DEALLOCATE (fixed) 1053 DEALLOCATE (rho_ij) 1054 DEALLOCATE (d_ij) 1055 DEALLOCATE (r_ij) 1056 DEALLOCATE (at_row) 1057 1058 END SUBROUTINE construct_initial_hess 1059 1060! ************************************************************************************************** 1061!> \brief ... 1062!> \param r1 ... 1063!> \param i ... 1064!> \param j ... 1065!> \param d ... 1066!> \param rho ... 1067!> \return ... 1068! ************************************************************************************************** 1069 FUNCTION dist_second_deriv(r1, i, j, d, rho) RESULT(deriv) 1070 REAL(KIND=dp), DIMENSION(3) :: r1 1071 INTEGER :: i, j 1072 REAL(KIND=dp) :: d, rho, deriv 1073 1074 deriv = 0.45_dp*rho*(r1(i)*r1(j))/d**2 1075 END FUNCTION 1076 1077! ************************************************************************************************** 1078!> \brief ... 1079!> \param r_ij ... 1080!> \param d_ij ... 1081!> \param rho_ij ... 1082!> \param idir ... 1083!> \param jdir ... 1084!> \param iat_der ... 1085!> \param jat_der ... 1086!> \param natom ... 1087!> \return ... 1088! ************************************************************************************************** 1089 FUNCTION angle_second_deriv(r_ij, d_ij, rho_ij, idir, jdir, iat_der, jat_der, natom) RESULT(deriv) 1090 REAL(KIND=dp), DIMENSION(:, :, :) :: r_ij 1091 REAL(KIND=dp), DIMENSION(:, :) :: d_ij, rho_ij 1092 INTEGER :: idir, jdir, iat_der, jat_der, natom 1093 REAL(KIND=dp) :: deriv 1094 1095 INTEGER :: i, iat, idr, j, jat, jdr 1096 REAL(KIND=dp) :: d12, d23, d31, D_mat(3, 2), denom1, & 1097 denom2, denom3, ka1, ka2, ka3, rho12, & 1098 rho23, rho31, rsst1, rsst2, rsst3 1099 REAL(KIND=dp), DIMENSION(3) :: r12, r23, r31 1100 1101 deriv = 0._dp 1102 IF (iat_der == jat_der) THEN 1103 DO i = 1, natom - 1 1104 IF (rho_ij(iat_der, i) .LT. 0.00001) CYCLE 1105 DO j = i + 1, natom 1106 IF (rho_ij(iat_der, j) .LT. 0.00001) CYCLE 1107 IF (i == iat_der .OR. j == iat_der) CYCLE 1108 IF (iat_der .LT. i .OR. iat_der .GT. j) THEN 1109 r12 = r_ij(iat_der, i, :); r23 = r_ij(i, j, :); r31 = r_ij(j, iat_der, :) 1110 d12 = d_ij(iat_der, i); d23 = d_ij(i, j); d31 = d_ij(j, iat_der) 1111 rho12 = rho_ij(iat_der, i); rho23 = rho_ij(i, j); rho31 = rho_ij(j, iat_der) 1112 ELSE 1113 r12 = r_ij(iat_der, j, :); r23 = r_ij(j, i, :); r31 = r_ij(i, iat_der, :) 1114 d12 = d_ij(iat_der, j); d23 = d_ij(j, i); d31 = d_ij(i, iat_der) 1115 rho12 = rho_ij(iat_der, j); rho23 = rho_ij(j, i); rho31 = rho_ij(i, iat_der) 1116 END IF 1117 ka1 = 0.15_dp*rho12*rho23; ka2 = 0.15_dp*rho23*rho31; ka3 = 0.15_dp*rho31*rho12 1118 rsst1 = DOT_PRODUCT(r12, r23); rsst2 = DOT_PRODUCT(r23, r31); rsst3 = DOT_PRODUCT(r31, r12) 1119 denom1 = 1.0_dp - rsst1**2/(d12**2*d23**2); denom2 = 1.0_dp - rsst2**2/(d23**2*d31**2) 1120 denom3 = 1.0_dp - rsst3**2/(d31**2*d12**2) 1121 denom1 = SIGN(1.0_dp, denom1)*MAX(ABS(denom1), 0.01_dp) 1122 denom2 = SIGN(1.0_dp, denom2)*MAX(ABS(denom2), 0.01_dp) 1123 denom3 = SIGN(1.0_dp, denom3)*MAX(ABS(denom3), 0.01_dp) 1124 D_mat(1, 1) = r23(idir)/(d12*d23) - rsst1*r12(idir)/(d12**3*d23) 1125 D_mat(1, 2) = r23(jdir)/(d12*d23) - rsst1*r12(jdir)/(d12**3*d23) 1126 D_mat(2, 1) = -r23(idir)/(d23*d31) + rsst2*r31(idir)/(d23*d31**3) 1127 D_mat(2, 2) = -r23(jdir)/(d23*d31) + rsst2*r31(jdir)/(d23*d31**3) 1128 D_mat(3, 1) = (r31(idir) - r12(idir))/(d31*d12) + rsst3*r31(idir)/(d31**3*d12) - & 1129 rsst3*r12(idir)/(d31*d12**3) 1130 D_mat(3, 2) = (r31(jdir) - r12(jdir))/(d31*d12) + rsst3*r31(jdir)/(d31**3*d12) - & 1131 rsst3*r12(jdir)/(d31*d12**3) 1132 IF (ABS(denom1) .LE. 0.011_dp) D_mat(1, 1) = 0.0_dp 1133 IF (ABS(denom2) .LE. 0.011_dp) D_mat(2, 1) = 0.0_dp 1134 IF (ABS(denom3) .LE. 0.011_dp) D_mat(3, 1) = 0.0_dp 1135 deriv = deriv + ka1*D_mat(1, 1)*D_mat(1, 2)/denom1 + & 1136 ka2*D_mat(2, 1)*D_mat(2, 2)/denom2 + & 1137 ka3*D_mat(3, 1)*D_mat(3, 2)/denom3 1138 1139 END DO 1140 END DO 1141 ELSE 1142 DO i = 1, natom 1143 IF (i == iat_der .OR. i == jat_der) CYCLE 1144 IF (jat_der .LT. iat_der) THEN 1145 iat = jat_der; jat = iat_der; idr = jdir; jdr = idir 1146 ELSE 1147 iat = iat_der; jat = jat_der; idr = idir; jdr = jdir 1148 END IF 1149 IF (jat .LT. i .OR. iat .GT. i) THEN 1150 r12 = r_ij(iat, jat, :); r23 = r_ij(jat, i, :); r31 = r_ij(i, iat, :) 1151 d12 = d_ij(iat, jat); d23 = d_ij(jat, i); d31 = d_ij(i, iat) 1152 rho12 = rho_ij(iat, jat); rho23 = rho_ij(jat, i); rho31 = rho_ij(i, iat) 1153 ELSE 1154 r12 = r_ij(iat, i, :); r23 = r_ij(i, jat, :); r31 = r_ij(jat, iat, :) 1155 d12 = d_ij(iat, i); d23 = d_ij(i, jat); d31 = d_ij(jat, iat) 1156 rho12 = rho_ij(iat, i); rho23 = rho_ij(i, jat); rho31 = rho_ij(jat, iat) 1157 END IF 1158 ka1 = 0.15_dp*rho12*rho23; ka2 = 0.15_dp*rho23*rho31; ka3 = 0.15_dp*rho31*rho12 1159 rsst1 = DOT_PRODUCT(r12, r23); rsst2 = DOT_PRODUCT(r23, r31); rsst3 = DOT_PRODUCT(r31, r12) 1160 denom1 = 1.0_dp - rsst1**2/(d12**2*d23**2); denom2 = 1.0_dp - rsst2**2/(d23**2*d31**2) 1161 denom3 = 1.0_dp - rsst3**2/(d31**2*d12**2) 1162 denom1 = SIGN(1.0_dp, denom1)*MAX(ABS(denom1), 0.01_dp) 1163 denom2 = SIGN(1.0_dp, denom2)*MAX(ABS(denom2), 0.01_dp) 1164 denom3 = SIGN(1.0_dp, denom3)*MAX(ABS(denom3), 0.01_dp) 1165 D_mat(1, 1) = r23(idr)/(d12*d23) - rsst1*r12(idr)/(d12**3*d23) 1166 D_mat(2, 1) = -r23(idr)/(d23*d31) + rsst2*r31(idr)/(d23*d31**3) 1167 D_mat(3, 1) = (r31(idr) - r12(idr))/(d31*d12) + rsst3*r31(idr)/(d31**3*d12) - & 1168 rsst3*r12(idr)/(d31*d12**3) 1169 IF (jat .LT. i .OR. iat .GT. i) THEN 1170 D_mat(1, 2) = (r12(jdr) - r23(jdr))/(d12*d23) + rsst1*r12(jdr)/(d12**3*d23) - & 1171 rsst1*r23(jdr)/(d12*d23**3) 1172 D_mat(2, 2) = r31(jdr)/(d23*d31) - rsst2*r23(jdr)/(d23**3*d31) 1173 D_mat(3, 2) = -r31(jdr)/(d31*d12) + rsst3*r12(jdr)/(d31*d12**3) 1174 ELSE 1175 D_mat(1, 2) = -r12(jdr)/(d12*d23) + rsst1*r23(jdr)/(d12*d23**3) 1176 D_mat(2, 2) = (r23(jdr) - r31(jdr))/(d23*d31) + rsst2*r23(jdr)/(d23**3*d31) - & 1177 rsst2*r31(jdr)/(d23*d31**3) 1178 D_mat(3, 2) = r12(jdr)/(d31*d12) - rsst3*r31(jdr)/(d31**3*d12) 1179 END IF 1180 IF (ABS(denom1) .LE. 0.011_dp) D_mat(1, 1) = 0.0_dp 1181 IF (ABS(denom2) .LE. 0.011_dp) D_mat(2, 1) = 0.0_dp 1182 IF (ABS(denom3) .LE. 0.011_dp) D_mat(3, 1) = 0.0_dp 1183 1184 deriv = deriv + ka1*D_mat(1, 1)*D_mat(1, 2)/denom1 + & 1185 ka2*D_mat(2, 1)*D_mat(2, 2)/denom2 + & 1186 ka3*D_mat(3, 1)*D_mat(3, 2)/denom3 1187 END DO 1188 END IF 1189 deriv = 0.25_dp*deriv 1190 1191 END FUNCTION angle_second_deriv 1192 1193END MODULE bfgs_optimizer 1194